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Abstract 


In  this  project,  computer  simulation  and  optical  modeling  of  laser  beam  propagation  through 
the  turbulent  atmosphere,  as  well  as  development  of  techniques  using  photorefractive  crystals 
to  mitigate  phase  distortions  in  laser  beams,  have  been  done.  In  the  framework  of  the  first 
direction  the  mathematical  methods  and  computational  schemes  based  on  split-step  operators, 
phase  screen  model  and  the  Monte  Carlo  method  have  been  elaborated.  Spatial  statistics  of 
light  field  in  laser  beam  has  been  studied  in  relation  to  inner  and  outer  scales  for  different 
models  of  atmospheric  turbulence.  Regimes  of  weak,  moderate  and  strong  fluctuations  have 
been  considered.  In  the  framework  of  the  second  direction  the  possibility  of  simulation  of 
double  pass  and  anisoplanatic  effects  by  means  of  few  phase  screens  has  been  studied.  An 
experimental  set-up  for  optical  modeling  anisoplanatic  effects  by  the  use  of  dynamic  phase 
modulator  has  been  designed.  A  method  of  generation  of  random  optical  field  with  variable 
correlation  function  has  been  proposed  and  tested.  In  the  framework  of  the  third  direction  the 
problem  of  mitigation  of  distorted  optical  signals  in  photorefractive  crystals  has  been  studied. 
An  optimal  effective  operating  range  of  one-way  system,  based  on  nonlinear  interaction  of 
distorted  signal  with  pumping  formed  by  spatial  filtering  of  the  signal ,  has  been  found.  Using 
the  criterion  of  maximal  mitigation  of  phase  and  amplitude  distortions,  the  schemes  of  two- 
and  four-beam  interaction  in  InP:Fe  have  been  optimized. 
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Part  I 

Numerical  Simulation  of  Laser  Beam  Propagation  Through 
Atmospheric  Turbulence 


1.  Introduction 

When  the  laser  beam  propagates  through  the  atmosphere  it  undergoes  distortions,  caused  by 
fluctuations  in  the  refractive  index  in  turbulence.  Random  blurring  of  the  beam  and  the  loss  of  the  spatial 
coherence  are  some  of  the  factors  which  limit  the  effectiveness  of  the  location,  probing  and  imaging  systems  in 
the  atmosphere. 

Lately  the  Monte  Carlo  method  has  been  intensively  developed  in  application  to  the  tlieoretical 
investigation  of  the  laser  beam  propagation  through  the  turbulent  atmosphere.  This  method  has  considerable 
advantages  over  the  analytical  methods,  which  commonly  use  one  or  another  set  of  approximations^  Rytov 
theory  is  applicable  in  the  regime  of  weak  fluctuations,  when  the  scintillation  index  in  the  plane  wave  satisfies 

the  condition  pi  «  1.  The  equations  for  the  field  moments  are  valid  without  restrictions  on  the  strength  of 
fluctuations^.  However,  for  the  kurtosis  (the  fourth  moment),  which  describes  irradiance  variance  and 
covariance,  only  approximate  solutions  may  be  obtained.  The  generalized  Huygens-Fresnel’s  approach  allows 

to  obtain  the  mean  irradiance,  but  the  solution  for  the  kurtosis  diverges  with  the  increasing  Rytov  variance  Pi . 
The  phase  approximation  in  Huygens-Kirchhoffs  method  gives  asymptotically  strict  solutions  for  the 

fluctuations  of  the  focused  beam  in  the  regions  Pl  «\  and  Pi  »\  and  for  the  fluctuations  of  the  collimated 
beam  only  in  the  near-field  region^.  The  interesting  results  have  been  obtained  from  the  theory  of  effective 
beam  parameters,  developed  in  Ref,  4.  The  Monte  Carlo  method  is,  in  principle,  free  from  any  limitations.  Its 
accuracy  is  defined  only  by  computing  resources. 

Optics  of  the  atmosphere  usually  assumes  two  approaches  in  the  development  of  the  Monte  Carlo 
method.  One  of  them  is  known  as  corpuscular  approach  and  the  other  as  wave  approach^*^  (Fig.  1.1a).  In 
corpuscular  approach  the  propagation  of  the  light  wave  through  the  atmosphere  is  considered  as  a  stochastic 
process  of  the  photon  scattering  by  the  molecules  of  the  air  components^  From  the  ensemble  of  several 
thousands  of  calculated  trajectories  the  angular  distribution  and  polarization  of  the  scattered  radiation®  as  well 
as  the  visibility  of  the  object  in  the  atmosphere^  and  the  quality  of  imaging  through  the  disperse  medium^°  may 
be  obtained. 

The  wave  approach  is  based  on  the  model  of  phase  screens  which  imitate  fluctuations  of  the 
permittivity  in  the  turbulent  atmosphere.  The  propagation  of  the  laser  radiation  is  considered  as  a  process  of 
successive  scattering  of  the  light  wave  by  the  phase  screens  (Fig.  1.1b).  The  numerical  solution  to  the  problem 
of  the  wave  scattering  yields  the  realization  of  the  light  field  in  the  receiver  plane.  Physically,  this  realization 
is  equivalent  to  the  distribution  of  the  light  field  registered  in  the  receiver  plane.  A  series  of  such  solutions 
obtained  for  different  sets  of  phase  screens  forms  an  ensemble  of  realizations  of  the  light  field.  Statistical 
characteristics  of  the  laser  beam  in  the  atmosphere  are  calculated  by  averaging  over  the  ensemble  of  the 
obtained  realizations. 

The  Monte  Carlo  method  allows  on  the  basis  of  a  unified  approach  to  investigate  various  statistical 
quantities,  associated  with  beam  propagation  in  turbulence,  as  well  as  various  models  of  atmospheric 
turbulence  and  conditions  of  propagation. 

The  Monte  Carlo  method  based  on  the  phase  screen  model  is  widely  used  for  the  treatment  of  the 
problems  of  the  wave  propagation  in  various  media.  The  method  is  applied  to  the  analysis  of  the  fluctuations 
which  are  observed  in  many  situations  of  scientific  interest,  such  as  electromagnetic  waves  through 
atmospheres,  ionospheres,  acoustic  waves  in  the  ocean,  seismic  waves  in  the  eartlFh 

The  present  paper  is  devoted  to  computer  simulations  of  the  statistics  of  the  laser  beam  in  the 
turbulent  atmosphere  on  the  basis  of  the  phase  screen  model.  This  paper  has  two  main  parts.  The  first  part 
(sections  2-4)  will  describe  methodological  aspects  of  the  Monte  Carlo  method  based  on  the  phase  screen 
model.  We  will  consider  the  justification  of  the  Monte  Carlo  metliod  and  the  phase  screen  model,  the  methods 
of  generating  the  phase  screens,  which  may  adequately  represent  atmospheric  turbulence,  and  methods  to  solve 
the  problem  of  scattering  of  the  radiation  by  a  set  of  phase  screens. 


6 


s  =  1 . S  ° 

Fig.  1.1 

The  Monte  Carlo  method  in  a  random  media,  a)  The  corpuscular  approach,  b)  The  wave  approach  on 
base  of  phase  screen  model. 


The  second  part  (sections  5-6)  will  discuss  the  results  of  the  Monte  Carlo  simulations.  We  will 
consider  the  effect  of  the  large-scale  turbulence  on  random  blur  and  wander  of  the  laser  beam  and  the  effect  of 
the  small-scale  turbulence  on  spatial  statistics  of  the  light  field  at  the  end  of  the  atmospheric  path.  The  effect  of 
the  model  of  atmospheric  turbulence  on  statistical  characteristics  of  the  beam  is  estimated.  The  results  obtained 
by  Monte  Carlo  simulations  and  analytical  methods  are  compared. 

2.  Justification  of  the  Monte  Carlo  method  based  on  the  phase  screen  model 

2.1,  Pambolic  equation 

The  wave  approach  in  the  Monte  Carlo  method  is  based  on  the  parabolic  equation  for  the  complex 
amplitude  of  the  light  field  E{p,z),  propagating  in  a  randomly  inhomogeneous  medium: 

lik.  —  =  ^,E  +  lI^^7i{p,z)E,  (2.1) 

dz 

where  k^-lni  wave  number,  n{p,z)  is  a  three-dimensional  field  of  the  refractive  index  fluctuations  in 
the  medium. 

Parabolic  equation  (2.1)  is  valid  under  the  following  assumptions^^: 

-  spatial  scales  of  fluctuations  in  the  random  medium  1^  and  in  the  light  field  are  much  larger  than  the 
wavelength  X  - 

4,4  »;i  (2.2) 

-  the  Fresnel  approximation  may  be  used: 

(2.3) 

-  the  energy  losses  due  to  the  backward  scattering  of  the  field  by  inhomogeneities  of  the  medium  are  small: 

2k 

n^lc^z  |0^(/c)A:(yA:«  1  (2.4) 

kyfl 

When  the  light  wave  propagates  through  the  atmosphere  the  conditions  (2.2)-(2.4)  are  always 

satisfied. 
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The  simullaneous  fulfillment  of  the  conditions  (2.2)-(2.4)  means  that  the  scattered  radiation 
propagates  in  the  forward  direction  and  the  backward  scattering  may  be  neglected.  This  fact  allows  to 

introduce  the  propagation  direction  z,  along  which  the  light  field  amplitude  variations  are  observed. 

The  evolution  of  the  radiation  along  the  coordinate  z  obeys  the  principle  of  dynamic  causality^^.  According  to 
it,  the  propagation  of  the  light  beam  is  Markov  process  of  successive  passage  of  the  slabs  of  randomly 
inhomogeneous  atmosphere. 

The  representation  of  continuum  as  a  set  of  slabs  has  been  used  in  the  local  method  of  small 
perturbations^^.  From  this  method  the  equations  for  the  field  moments  may  be  immediately  obtained. 

2,2.  Model  of  phase  screens 

The  further  development  of  the  stratified  representation  of  continuum  is  a  model  of  phase  screens.  In 
this  model  each  slab  of  length  Az  is  effectively  compressed  to  a  phase  screen  located  at  the  slab  center.  At  a 
distance  Az/2  to  the  right  and  to  the  left  from  the  phase  screen  only  diffraction  in  free  space  is  considered 
(Fig.  1.1b).  Randomly  inhomogeneous  continuum  is  replaced  by  a  set  of  phase  screens  located  along  the 
propagation  direction.  The  distance  between  the  screens  is  Az,  the  first  and  the  last  screen  arc  located  at  a 
distance  Az/  2  from  the  transmitter  and  receiver  planes  respectively. 

Thus,  fluctuations  in  the  refractive  index  of  the  atmosphere  induce  fluctuations  in  the  phase  of  the 

light  field  (p).  These  fluctuations  are  located  on  the  infinitely  thin  phase  screens  in  the  planes  4+1/2. 

Between  the  screens  phase  fluctuations  transform  to  amplitude  fluctuations  owing  to  diffraction. 

Such  a  set  of  a  finite  number  of  phase  screen  is  a  good  model  of  randomly  inhomogeneous  continuum, 

if  the  length  Az  between  the  screens  is  small  in  comparison  with  the  diffraction  length  and  extinction 
length  a~^ . 


Az<min{Z</,a”^},  Az=4‘“4-i» 

where  L^^krl  -  is  diffraction  length  corresponding  to  the  characteristic  scale  of  the  light  field 
inhomogeneity,  a  is  the  extinction  coefficient. 

When  the  condition  (2.5)  is  satisfied  the  phase  fluctuation  on  the  phase  screen  can  be  obtained  from 
the  geometrical  optics  approximation: 

=  ^  (2.6) 

If  the  length  of  a  slab  Az  exceeds  the  outer  scale  of  turbulence 

Az»L^,  (2.7) 

then  the  model  of  J -correlated  relative  to  z  fluctuations  of  the  refractive  index  is  valid.  For  statistically 
isotropic  field  of  atmospheric  turbulence  the  spatial  spectrum  of  the  random  phase  g){p)  on  the  screen 

may  be  expressed  as  follows: 

F^{kJ  =  27±^Az<l>,{kj_,0),  (2.8) 

where  is  a  three  -  dimensioiral  spectrum  of  the  refractive  index  fluctuations  in  the  atmosphere. 

The  variance  aj‘  of  the  phase  fluctuations  on  the  screen  is  given  by; 

CO 

=  2  ,0)x:£fc  (2.9) 

0 

For  Gaussian  field  h{p,z)  the  condition  (2.7)  means  that  phase  fluctuations  on  different  screens 
are  statistically  independent: 

=  (2.10) 

When  considering  the  wind  flow,  one  commonly  uses  the  hypothesis  of  the  "frozen  turbulence".  In 
this  case  the  field  (p{p)  satisfies  the  condition: 
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(2.11) 


where  Vi  is  the  wind  flow  velocity  in  the  plane  of  the  screen. 

2J.  The  origin  of  the  phase  screen  model 

The  simplest  model  containing  one  phase  screen  reflects  the  main  properties  of  the  wave  in  randomly 
inhomogeneous  medium.  This  model  allows  to  qualitatively  estimate  statistical  properties  of  the  wave  in 
random  medium.  Presumably,  the  model  of  phase  screen  was  first  used  in  Ref.  14  for  the  analysis  of  the  wave 
diffraction  in  a  thick  slab  of  randomly  inhomogeneous  medium.  The  theoiy  of  the  wave  diffraction  behind  a 
thin  screen  with  random  phase  and  amplitude  fluctuations  was  thoroughly  considered  in  Ref  15.  The  physical 
situation,  which  is  most  adequately  described  by  one  phase  screen,  is  scattering  of  the  wave  in  a  thin  slab, 
when  fluctuations  of  the  light  field  are  formed  behind  the  slab.  In  particular,  this  situation  may  occur  when  the 
wave  propagates  through  the  ionosphere^^  and  interplanetary  plasma^^. 

The  error  of  the  representation  of  the  continuum  by  one  phase  screen  has  been  studied  in  Ref  17  on 
the  basis  of  numerical  solution  of  the  for  the  field  kourtosis.  When  the  length  of  the  medium  exceeds  the  outer 

scale  Lq,  one  phase  screens  significantly  overestimates  irradiance  fluctuations  in  the  absence  of  saturation. 

The  error  of  the  model  essentially  decreases  if  the  screen  is  located  at  the  center  of  the  path. 

The  first  statistical  trials  with  the  waves  scattered  by  the  phase  screen  have  been  done  in  Ref  18.  The 
authors^*  discuss  methodological  questions  of  spectral  analysis  of  stochastic  waves  diffraction  and  illustrate  the 
methods  by  the  example  of  propagation  through  a  thin  layer  of  anisotropic  ionosphere.  The  distribution 
function  of  the  irradiance  fluctuations  in  the  plane  wave  behind  a  thin  turbulent  layer  was  studied  in  Ref  19  on 
the  basis  of  the  model  with  one  phase  screen.  Ten  thousand  realizations  showed  that  Rice-Nakagami  and  log¬ 
normal  distributions  may  not  be  applied  to  the  description  of  irradiance  fluctuations. 

The  model  based  on  a  set  of  phase  screens  is  substantially  more  reliable  for  the  reproduction  of 
randomly  inhomogeneous  continuum.  The  application  of  this  model  is  equivalent  to  the  representation  of  the 
kurtosis  of  the  field  as  a  convolution  for  the  continuum^^.  In  the  plane  wave  coherence  function  integration  in 
the  exponent  is  replaced  by  summation  when  the  phase  screen  method  is  applied^^ 

The  application  of  the  phase  screen  model  to  analytical  study  of  the  wave  propagation  in  randomly 
inhomogeneous  media  is  thoroughly  discussed  in  Refs.  6, 11. 

Statistical  study  of  the  wave  propagation  in  randomly  inhomogeneous  media  on  the  basis  of  the  phase 
screen  model  implies  the  application  of  the  split-step  method  to  solve  stochastic  parabolic  equation  (2.1). 

2,4.  Phase  screen  inodei  and spUt-step  method 

In  fact,  the  phase  screen  model  for  randomly  inhomogeneous  medium  is  physical  interpretation  of 
split-step  method,  which  is  well  known  in  computational  mathematics^^.  In  this  method  the  differential 
equation  with  some  operator,  which  can  be  written  as  a  linear  sum  of  several  pieces,  at  each  step  Az  is 
replaced  by  a  sequence  of  equations  corresponding  to  the  pieces  in  the  sum.  This  approach  allows  to  use  the 
most  effective  differencing  schemes  for  each  equation  from  the  sequence.  It  has  been  shown^^  that  the  split- 
step  method  allows  to  reduce  computing  time  by  as  much  as  30  percents  in  comparison  with  the  time  needed 
for  the  direct  solution  of  the  stochastic  equation  (2.1)  for  randomly  inhomogeneous  medium.  When  handling 
diffraction  problems  Fast  Fourier  Transform  algorithm  is  found  to  especially  increase  effectiveness  of  the  split- 
step  method  in  the  conditions  of  strong  fluctuations. 

The  idea  to  use  the  split-step  method  for  treatment  of  the  wave  propagation  in  randomly 
inhomogeneous  and  nonlinear  media  was  first  suggested  in  Ref  24.  One  of  the  first  practical  applications  of 
this  step  method  was  associated  with  investigation  of  the  effect  of  the  internal  waves  on  distortions  of  acoustic 
signal  in  undenvaler  sound  channeF^.  In  optics  of  the  turbulent  atmosphere  the  split-step  method  has  become 
widespread  after  Ref  26  had  been  published.  In  Ref  26  separate  realizations  of  a  powerful  laser  beam 
traveling  through  the  atmosphere  were  obtained. 

The  current  state  concept  of  the  split-step  method  and  ways  for  constructing  the  splitted  operators  for 
the  problems  of  the  optical  wave  propagation  in  randomly  inhomogeneous  medium  are  reviewed  in  Ref  27. 

For  the  parabolic  equation  (2.1)  the  split-step  method  applied  to  a  slab  Az  of  the  medium  yields  the 
following  set  of  equations: 


=  ze[^,z^yi],  fox  E^{p,z,)  =  E{p,z,)\ 

oz 


(2.12) 
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(2.13) 


lika-^=21^^hE^,  2-e[z„z^,]  for  eXp,z,)  =  E^{p,z^^^,^)\ 
oz 

2tk^^-  =  l^iEj,  zB\z^^n,Zs*^,  for  i'^U-^^i/z)  =  ^f(A^i+i)-  (214) 

The  solution  of  the  last  equation  is  the  final  solution  of  the  equation  (2.1)  in  the  plane  ZsW. 
EAp>Zs^i)  =  E{p,z,^X 

Equations  (2.12)  and  (2.14)  describe  the  diffraction  of  the  wave  in  homogeneous  field  at  a  distance 
Az/2  to  the  left  and  to  the  right  of  the  screen,  located  in  the  plane  Equation  (2.13)  describes  the 
scattering  of  the  wave  by  inhomogeneities  in  geometric  optics  approximation.  Its  solution  is 

=  (2,15) 

where  ^^+1/2  (p)  is  the  phase  perturbation  on  the  screen  (2,6). 

In  the  next  slab  -^+2]  we  sequentially  solve  equations  describing  diffraction  and  scattering  of  the 
wave  by  a  phase  screen.  This  procedure  repeats  over  and  over  again  from  transmitter  plane  2=0  to  receiver 
plane  z=z^. 

It  is  common  practice  to  combine  equations  describing  diffraction  in  two  adjacent  slabs.  As  a  result 

diffraction  of  the  wave  is  considered  on  the  interval  [-^^_i/2 +0,^5+i/2  i-^-  between  two  neighboring 

screens. 


Z5.  Statistical  trials  in  optics  of  the  turbulent  atmosphere 

First  papers  dealing  with  the  Monte  Carlo  method  based  on  the  phase  screen  approach  are  related  to 
investigations  of  statistical  characteristics  of  the  powerful  laser  beams  in  the  turbulent  atmosphere.  In  Ref.  28  a 
time-dependent  calculation  of  spatial  statistics  of  irradiance  fluctuations  of  the  beam  propagating  through  the 
atmosphere  under  the  conditions  of  thermal  self-action  has  been  done.  Irradiance  fluctuations  of  the  laser  beam 
in  the  presence  of  wind  are  considered  in  Ref.  29.  The  effect  of  the  wind  velocity  fluctuations  on  random 
motion  and  blurring  of  the  beam  propagating  in  the  presence  of  nonlinear  refraction  is  discussed  in  Refs. 
30,31. 

Recently  much  attention  has  been  paid  to  the  linear  problems  of  statistical  optics  of  the  turbulent 
atmosphere.  Here  the  propagation  is  described  by  the  equation  (2.1).  A  series  of  Refs.  32-34  have  been  devoted 
to  the  investigation  of  the  plane  and  spherical  light  waves  statistics.  In  Ref  32  strong  fluctuations  of  the  plane 
wave  irradiance  are  considered,  the  Ref  33  investigates  the  convergence  of  asymptotic  solutions  for  irradiance 

moments  in  the  regimes  of  strong  fil  and  complete  pl  ^  10  saturation  of  fluctuations.  In  Ref  34  simple 
empirical  formulas  have  been  obtained  for  the  dependence  of  irradiance  variance  on  the  inner  scale  in  the 
regime  of  strong  focusing.  The  Ref  35  is  devoted  to  investigation  of  the  distribution  function  of  irradiance 
fluctuations  in  the  plane  and  spherical  waves. 

3.  Computer  simulations  of  the  phase  screens 

3.1.  Problem  of  scales  in  atmospheric  turbulence 

According  to  the  modern  concepts  the  microstructure  of  atmospheric  turbulence  is  determined  by  the 
cascade  process  of  fragmentation  of  eddies  that  arrise  in  the  wind  flow^^.  The  size  of  the  large-scale  eddies 

/q  «r«L^  is  comparable  to  the  characteristic  scale  of  the  entire  wind  flow.  The  size  of  the  small-scale 
eddies  4  corresponds  to  the  size  of  perturbations  at  which  the  dissipation  of  energy  through  viscous  effects 
becomes  important.  The  inertial  subrange  of  atmospheric  turbulence  is  bounded  by  the  outer  and  inner  /q 
scales 


l^«r«I^  (3.1) 

Structure  functions  of  the  random  temperature  field  and  consequently  the  refractive  index  field  obey  within  the 
inertial  subrange  the  two-third  law^: 
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Io«r«I„  (3-2) 

where  is  the  structure  constant  of  the  refractive  index  fluctuations  in  the  atmosphere.  This  structure 

function  corresponds  to  the  three-dimensional  Kolmogorov  spectrum  Of  (*:)  that  has  a  power  dependence  on 
2 

the  wave  number  ic  : 

0/(*:)  =  0.033C/k""'\  k,  «k«  k,  (3-3) 

where  =  27i  /  /„ ,  Ko  =  27t  /  4  are  the  boundaries  of  the  spatial  spectrum  corresponding  to  the  inertial 

subrange.  ^  j 

The  Kolmogorov  spectrum  satisfactorily  describes  the  fluctuations  of  the  refractive  index  in  the 

atmosphere  within  the  inertial  subrange.  However,  the  behaivor  of  the  fluctuations  at  the  boundaries  and 
outside  of  the  inertial  subrange  requires  some  additional  assumptions.  The  following  models  take  into  account 
the  effect  of  the  boundaries  in  an  explicit  form: 
the  Tatarskii  spectrum  considers  the  effect  of  the  inner  scale 


the  von  Karman  spectrum  considers  the  effect  of  the  outer  scale 

<l)/^(K)  =  0.033q,'(x'  +Aro')-”'^ 

the  modified  von  Karman 

o/"^(k)  =  0.033  C/ (x:^ 

and  the  Andrews  spectra 


0/(k)  =  0.033G'  1  +  1.8021 


-0.2541  — 


r.  =  5.92/.4 

(3.4) 

(3.5) 

a:„  =  5.92/^ 

(3.6) 

- - - ,  =  3.3/^ 

(/e  +  ^o')"'‘ 

(3.7) 

consider  the  effect  of  both  inner  and  outer  scales 

Typically,  the  scale  sizes  of  turbulent  eddies  vary  from  an  outer  scale  of  tens  of  meters  down  to  an 
inner  scale  of  just  a  few  millimeters.  In  other  words,  the  range  of  spatial  perturbations  of  the  refractive  index 
in  the  atmosphere  is  extremely  wide:  the  ratio  between  the  iimer  and  the  outer  scales  may  reach  several  orders 

of  magnitude.  In  numerical  implementations,  K;  may  be  often  ignored,  because  the  ratio  of 

<l)j,(/Co,z)/Oj,(x:/,z)  or  the  dynamic  range  of  the  spectrum  to  be  simulated  is  several  orders  of  magnitude  as 

well.  Hence,  the  effect  of  the  roll-off  due  to  the  presence  of  in  the  model  will  most  likely  not  be  noticed. 

In  the  phase  screen  model  we  use  the  2D  power  spectrum  of  a  thin  slab  of  turbulence  in  a  plane 
orthogonal  to  the  direction  of  propagation  (the  zaxis).  This  spectrum  is  given  by  (3.8) 

/'(Ki)  =  2;rio^A20„(Ki,o)  (3.8) 

where  O„(ki,0)  is  2D  power  spectrum  of  the  refractive  index  fluctuations,  which  corresponds  to  one  of  the 
models  of  atmospheric  turbulence. 

3.2.  FFT-based phase  screen  generation 

The  traditional  method  of  generating  random  fields  with  the  known  autocorrelation  function  (the 
phase  screens  under  discussion  also  belong  to  this  case)  is  filtering  a  white  noise  process  with  square  root  of 
the  spectrum,  followed  by  an  inverse  Fourier  transform.  This  spectral  method  is  highly  efficient,  because  it 
allows  to  take  advantage  of  the  fast  Fourier  transform  (FFT)  algorithm. 

In  this  method  the  complex  field  is  represented  on  the  screen  as  a  sum  of  Fourier  harmonics 

with  random  complex  coefficients  (discrete  white  noise): 
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(3.9) 


=  ^F^{p,q)AK^AKy ,  =  expji-^j , 

where  ^pg^r/j^g  are  statistically  independent  random  numbers  which  are  either  Gaussian  with  parameters 
(0,1)  or  uniformly  distributed  on  the  interval  [-^6 /2,V6 /2].  In  the  last  case  the  generated  phase  screens 
have  Gaussian  statistics  owing  to  the  central  limit  theorem^^.  Through  F^{p,q)  =  a:^)  we  denote  the 

values  of  the  two-dimensional  spectrum  of  the  phase,  where  spatial  frequencies  a:^  and  Kg  are  defined  by 

P=  -NJ2,...NJ2-l, 

Kg  =  AKpj  q=-Nyl2,.„Nyl2-\,  ,  (3.10) 

Frequency  intervals  between  harmonics  Ak^  and  AKy  are  given  by 

AKx-2nl  A^y  AKy~2nl  Ay  (3.11) 

where  A^=  Nj,Ax  and  Ay=  NyAy  are  the  sizes  of  the  simulation  aperture.  Real  and  imaginary  parts  of  the 

generated  field  g){nyn^  form  two  statistically  independent  realizations  of  the  phase  screen 

{(p'(n,m)  =  Rc^n,w),  n=l,...,N^,  i2i=l,...,Ny] , 

[g)"{n,in)  =  lm<p{n,n^,  n  =  l,...,K„  (3.13) 

Algorithm  (3.9)  assumes  that  p  and  q  may  be  equal  to  zero.  In  this  case  the  weighting  coefficient  a^^is  equal 

to  the  value  of  zero  harmonic  of  the  spectrum.  For  the  Kolmogorov  (3.3)  and  Tatarskii  (3.4)  spectra  of 

atmospheric  turbulence  AJ,  co  as  a:  ^  0 .  Taking  into  accont  that  the  average  phase  delay  has  no  effect  on 

the  image  formation  process  or  a  beam  propagation,  the  origin  of  the  filter  function  can  be  set  to  zero,  =  0 . 

The  two-dimensional  field  of  the  phase  fluctuations  and  its  covariance  and  structure  functions  are 
represented  by  their  values  at  the  discrete  spatial  domain  poins  n,  m.  The  spectrum  of  the  field  is  represented 

by  its  values  at  the  discrete  spectral  domain  points  Kp^Kg.  Therefore,  the  phase  covariance  function 

and  its  spectrum  F^p,  g)  are  related  by  the  discrete  Fourier  transform 

B,{n,w)=  ""T  ""e”'  F{p,q)^KMy%.^%^  (3.14) 

Eq.  (3.14)  defines  the  simulated  (which  can  be  simulated  on  the  numerical  grid)  phase  covariance  function  for 
the  chosen  spectrum  of  atmospheric  turbulence  Fg,[K^yK^,  The  simulated  structure  function  is  given  by 

D^{n,in)  =  l[B,{0,Q)-B^{n,w))  (3.15) 

A  deviation  of  the  simulated  statistical  characteristics  from  the  desired  ones  forms  a  systematic  error 
of  the  phase  screens.  Typically,  phase  screens  are  evaluated  by  how  well  they  reproduce  the  desired  phase 
structure  function  for  the  given  turbulence  model.  For  the  Kolmogorov  turbulence  the  phase  structure  function 
has  the  form  (3.2).  When  discussing  methods  of  generating  the  phase  screens,  in  parallel  with  the  Kolmogorov 
spectrum  (3,3)  we  will  consider  the  von  Karman  spectrum  (3.5).  For  the  case  of  a  finite  outer  scale  the  phase 
structure  function  is  given  by 


where  Yy.^{rK^  is  Ihc  Neumann  function  of  the  order  of  1/3  of  the  imaginary  argument. 


(3.16) 
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Generation  of  random  screens  according  to  the  algorithm  (3.9)  can  be  easily  performed  using  the  FFT 
method.  However,  tins  fast  and  effective  approach  has  its  own  limitations.  The  minimum  and  ma.\imum 

spatial  frequencies  of  the  screens  generated  by  the  FFT  are  =  Ak  =  27:  /  >1  and  =  Mk  /  2  =  t:  /  Ax. 

Correspondingly,  the  minimum  representable  scale  size  is  =2Ax=27r/K^^  .  Thus,  if  one  wants  to 
reproduce  on  the  grid  the  minimum  inhomogeneity  of  the  order  of  10mm  and  if  the  number  of  grid  points  for 
one  of  the  transverse  coordinates  is  128,  than  the  maximum  reproducable  scale  on  the  grid  is  0.64  m.  In  most 

cases  this  is  not  enough  for  adequate  reproduction  of  the  large  scale  fluctuations  1^)  of  the  refractive  index 
since  typically  is  much  large  than  0.6m. 


Fig.  3.1.  Structure  function  of  the  phase 
fluctuations.  The  structure  constant 


=  10"^^cm■^^^  the  length  of  the  turbulent  slab 
Az-  1.06p  m,  the  minimum  representable  size 

7-0=0,1/27,  the  number  of  the  sample  points 
N-12d,  the  number  of  realizations  /W=100.  Monte 
Carlo  method  —  dotted  curve,  confidence  limits  — 
long  dashes,  analytical  dependence  -  dashed 
curve  (short),  simulated  (expected)  structure 
functions  —  solid  curve,  a  —  The  Kolmogorov 
spectrum,  b,  c  -  the  von  Karman  spectrum  with 

=50/77  (b)  and  4  m  (c). 

Fig.  3.1  shows  the  structure  functions  of 
phase  fluctuations  on  the  screen  obtained  by  the 
FFT-based  method.  The  averaging  has  been  done 


over  100  screens.  It  is  clearly  seen  that  analytical  dependences  lie  essentially  upper  than  the  curves  obtained 
through  Monte  Carlo  simulations  as  well  as  than  the  values  of  structure  function  expected  on  the  grid.  This  is 
valid  for  both  the  Kolmogorov  (Fig.  3.1a)  and  the  von  Karman  (with  large  outer  scale)  (Fig.  3.1b)  spectra.  A 
satisfactory  agreement  between  the  analytical  and  simulated  structure  fuinctions  is  obtained  only  when  the  step 
of  the  grid  becomes  comparable  with  the  outer  scale  of  turbulence  (Fig.  3.1c).  The  reason  for  the  discrepances 
shown  in  Figs.3.1  a,  and  b  is  that  the  tilt  and  other  manifestations  of  the  large-scale  turbulence  are  not  correct 
on  this  screen.  In  other  words,  the  grid  does  not  adequately  reproduce  the  low-frequency  characteristics  of 


atmospheric  turbulence. 
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J.J.  Subharmonks  method 

The  discussed  limitations  essentially  deteriorate  the  performance  of  the  FFT-based  methods  of 
generating  the  phase  screens.  In  a  series  of  papers^  some  attempts  have  been  done  to  modify  the  original 
method  so  that  it  could  take  into  account  the  effect  of  large-sc^e  phase  fluctuations  within  the  limited 
aperture. 

Following"^^  let  us  consider  the  possibilty  to  add  low-frequency  information  using  subharmonics.  The 

method  is  based  on  the  following  idea.  The  power  spectral  density  of  the  phase  fluctuations  in  the 

vicinity  of  zero  is  sampled  with  the  smaller  sampling  interval  than  at  the  rest  (high-frequency)  part  of  the 
spectrum.  Subsequent  application  of  this  procedure  leads  to  the  algorithm  of  increasing  (in  the  vicinity  of  zero) 
sampling  density  or  decreasing  sampling  interval  of  spatial  frequencies.  In  this  case  the  generation  of  the 
phase  screen  consists  of  several  steps. 

The  first  step  assumes  the  generation  of  the  phase  screen  (pj^(/7,zz7)  using  conventional  FFT-based 

method.  However,  before  making  inverse  Fourier  transformation  (3.9)  harmonics  F^{p,q)  are  modified  so 
that  the  points  at  indices  (±1,  0)  and  (0,  ±1)  are  scaled  by  a  factor  of  1/2  and  the  points  at  indices  (±1,  ±1)  are 
scaled  by  a  factor  of  3/4.  The  obtained  field  represents  the  high-frequency  part  of  the  phase 

fluctuations  on  the  screen. 

In  the  second  step  the  low-frequency  part  of  the  spectrum  is  formed  according  to  algorithm: 


''(p'  +  0.5)n]  ^  (q'  +  0.5)n'^ 


N, 


J 


(3.17) 


where  =3  +0.5)Aic^  and  =  3~^  (g' +0.5)Ak^.  The  values  of  harmonics  with  indices 

(-1,-1),  (0,-1),  (-1,0)  and  (0,0)  are  set  to  zero  for  all  y  except  for  J=  Nj.  The  number  Nj  defines  the  number 
of  iterations  of  increasing  sampling  density. 

For  each  iteration  sampling  interval  in  the  spatial  frequency  domain  decreases  by  a  factor  of  3.  The  larger  the 
number  of  iterations  the  closer  to  zero  frequency  are  harmonics  in  the  sum  (3.17). 

The  resulting  phase  screen  may  be  obtatined  as  a  sum  of  its  low-  and  high-frequency 

components,  cp(/?,/i^  =  (p^(/7,/z7)  +  9/(/7,/z^  . 

Fig.  3.2  demonstrates  the  results  of  statistical  processing  of  the  phase  screens  obtained  by 
subharmonics  method.  The  number  of  points  on  the  screen  is  the  same  as  in  the  FFT-based  method.  One  can 
see  that  the  increased  sampling  density  has  resulted  in  the  structure  functions  which  are  very  closed  to  those 
predicted  by  the  theory.  An  excellent  agreement  is  obtained  for  the  Karman  spectrum  with  the  outer  scale 

=  50/27. 

In  the  case  of  the  Kolmogorov  spectrum  the  results  of  Monte  Carlo  simulations  are  essentially  closer 
to  the  theoretically  predicted  as  compared  to  those  obtained  by  the  conventionsal  FFT  method  (Fig.  3.1a).  The 

agreement  may  be  further  improved  by  increasing  the  number  of  iterations  Nj.  However,  we  have  to  keep  in 


mind  that  with  increasing  number  of  iterations  the  round-off  error  also  increases  since  >  oo  when 

K-^0. 
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Fig.  3.2.  Structure  function  of  the  phase  fluctuations  on  the  screens  generated  by  subharmonic  metod.  Nj-A. 
Parameters  are  the  same  as  in  Fig.  3,1.  Solid  curve  —  analytical  dependence,  a  —  the  Kolmogorov  spectrum,  b 
-  the  von  Karman  spectrum  with  =  5Qm. 

A  limitation  of  this  method  is  its  low  efficiency.  Since  the  FFT  algorithm  is  not  used  when  the  low- 
frequency  part  of  the  spectrum  (3.17)  is  constructed,  the  computing  time  may  increase  by  an  order  of 
magnitude  in  comparison  with  the  FFT-based  method,  even  if  the  number  of  realization  is  not  too  large. 

3,4.  Modal  representation 

Along  with  the  spectral  method,  there  exist  alternative  techniques  for  generating  phase  screens.  One 
of  them  is  modal  representation.  In  contrast  to  the  spectral  method,  in  which  sines  and  cosines  form 
orthonormal  basis  for  the  decomposition  of  phase  fluctuations,  the  modal  method  employs  another  set  of 
orthonormal  functions.  For  phase  fluctuations  in  the  turbulent  atmosphere  the  most  appropriate  basis  is  formed 
by  Zemike  polynomials.  The  analytical  definition  of  these  polynomials  is  not  unique  since  there  is  a  choice 
with  respect  to  the  numbering  sequence  and  the  selection  of  an  appropriate  amplitude  factor.  We  will  use  the 
definition  convinient  from  the  point  of  view  of  a  statistical  analysis  ; 


y  =  ^/^TO^/2cos(nfi)l 

j  =  i?'(r)  V2  sHnB)  J  ,  (3.18) 

Zj  =  y[^]^(r),  m=0 

where 


_ (-l)"(n-5)! _ 

S  s![(n+/2j)/2-5]![(n-/2:^/2-5]! 


(3.19) 


The  values  of  n  and  m  are  always  integer  and  satisfy  /?,  /j-j/rj  =  even.  The  index  y  is  a  mode  ordering 
number  and  is  a  function  of  n  and  m. 

The  main  advantage  of  Zernike  polynomials  (3.18)  is  their  ability  to  represent  phase  fluctuations  in 
terms  of  classical  aberrations  observed  in  the  turbulent  atmosphere.  The  first  Zemike  polynomial 
reproduces  the  average  random  phase  growth  on  the  aperture  the  second  and  the  third  {Zi  and  describe 
random  tilt,  the  fourth  (^)  -  random  defocus.  The  Zernike  polynomial  expansion  of  random  phase  fluctuations 
on  the  screen  over  a  circle  of  radius  R  is  given  by: 

9(^,0)  =  Z5y^y(p;e)  (3.20) 

j 
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with  p  =  r/R  and  the  coefficients  aj  being  given  by 

aj  =  J  r/V (Hp)cp(%e)  Zj{p,%)  ,  (3.21) 

where 


p^i 

p>i 


(3,22) 


Statistical  properties  of  coefficients  Bj  depend  on  the  characteristics  of  the  wavefront  (p(i^,0)  to  be 

expanded.  In  the  case  under  discussion  aj  are  Gaussian  random  numbers  with  zero  mean  and  given  variance. 
The  covariance  matrix  of  these  coefficients  is  defined  by: 

Cj^  =<a-a^>  =  l  dpi  dp'  Mp)  <  q>(^.0)  cp(i?iP',0')  >  ^  (p'.G') .  (3.23) 

In  frequency  domain  we  have: 

C,y=<a;  >  =  ldi^dk'(^(Jr,ms(^/J^H^-^)QA^’'^'')’  (3.24) 

where  is  the  Fourier  transform  of  Zj{p,6).  The  function  <?y(ir,<j))  can  be  written  from  Eq.(3.18) 


Qodd 


j{m= 


Vrr  +  l 


QMA)  = 


nk 


'(-l)''^">''/”V2cos(/n|)) 

-  (_!)('■-”)« /”>7Jsin(/4)  , 
(-1)"'\  (m=0) 


(3.25) 


where  /[(a:)  is  the  /th  order  Bessel  fimction  of  the  first  kind. 

The  expressions  for  the  covariance  matrix  Cmay  be  obtained  analytically  only  for  the  Kolmogorov 
model  of  atmospheric  turbulence  (3.3)"^^  For  other  models  numerical  integration  of  the  expression  (3.24) 
should  be  performed. 

Although  the  Zernike  polynomials  form  an  orthonormal  basis,  it  follows  from  the  Eq.(3.23)  that  in  the 
general  case  coefficients  aj  are  not  statistically  independent.  The  latter  makes  the  generation  of  phase  screens 
using  Zernike  expansion  significantly  more  complicated.  The  requirement  of  statistical  independency  imposed 
on  coefficients  aj  in  combination  with  the  requirement  on  polynomials  to  be  orthonormal,  leads  to  the  system 
of  Karhunen-Loeve  equations: 


I  j  dp\  4>'2,(P,e)z,  (p',eo  =  5,^  ^ 

|Jc/pJ4D'z/p,0)<cp(;^,0)cp(i^^0O>z^  ^  ^ 

Thus,  the  Zernike  polynomials  are  not  orthonormal,  while  the  Karhunen-Loeve  functions  Zj  are  not 

convenient  since  they  cannot  be  expressed  analytically. 

In  the  present  paper  we  suggest  the  following  method  of  generating  random  fields  with  given  spatial 
spectrum.  The  method  is  based  on  the  Zernike  polynomial  expansion  of  the  Karhunen-Loeve  functions"^^.  Let 


^2 


\^Jj 


be  the  desired  vector  the  components  of  which  are  the  random  coefficients  for  the  Zernike 


polynomials.  The  covariance  matrix  C=  may  be  obtained  by  numerical  integration  of  the  Eq.(3.24) 
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the  Karhunen-Loeve  functions,  one  should  diagonalize  the  matrix  C.  It  is  known,  that  the  symmetric  quadratic 
matrix  may  be  diagonalized  by  similarity  transformation  S=U-C-lF,  where  is  a  unitary  matrix.  Here 
diagonal  elements  of  the  matrix  S  are  eigenvalues  of  the  matrix  Q  while  the  columns  of  the  matrix  U  are 
eigenvectors  of  the  initial  matrix.  However,  if  there  exists  indentical  eigenvalues  of  the  matrix  C,  not  any 
matrix  U  may  be  used  for  the  diagonalization  of  the  matrix  C,  In  our  case  the  modes  with  the  same  radial 

degrees  have  indentical  eigenvalues.  At  the  same  time  it  can  be  easily  verified  that  <  ay  ay  >  9^  0  if  the 

corresponding  polynomials  have  identical  angular  dependence.  This  fact  allows  to  group  the  coelficients, 
which  are  correlated,  in  clusters  and  apply  the  diagonalization  algorithm  within  clusters.  The  problem  of 
eigenvectors  disappears  automatically  since  the  polynomials  with  identical  angular  dependence  have  different 
radial  degrees. 

After  the  diagonalization  of  matrix  C  we  easily  obtain  the  vector  of  random  Karhunen-Loeve 
coefficients  B,  components  of  which  are  Gaussian  random  variables  with  zero  mean  and  variance,  gwen  by  the 
diagonal  matrix  S.  The  components  of  the  vector  A  can  be  computed  from  the  relation  A-U  -B.  These 

components  Bj  multiplied  by  {IR! where  75  =  1.68(CJPAz)-^''  is  the  Fried  radius,  give  rise  to  the 
weights  for  evaluating  the  Zernike  polynomials  with  correct  variance  and  covariance. 

Statistical  properties  of  the  screens,  obtained  from  the  modal  method  (3.20),  have  been  verified  by 
analysis  of  structure  functions  of  phase  fluctuations  on  these  screens.  Parameters  of  turbulence  and  grid  are  as 
follows: 

cm-2/3,  2=0.53  /mi,  A^50  m,  A^  =  A/^  =  128,  the  number  ofithe  Zernike  polynomials  /varied 

from  2  to  104.  The  averaging  has  been  performed  over  M=300  screens. 

In  order  to  make  sure  that  the  phase  screens  are  statistically  homogeneous  and  isotropic,  we  computed 
the  structure  function  in  several  cross-sections  of  the  screen.  Fig.  3.3  shows  three  structure  functions,  two  of 
which  have  been  computed  along  the  x-  and  7-  axes,  while  the  third  -  along  the  chord  7=1/3  i?).  The 
comparison  of  these  dependences  allows  to  conclude  that  statistical  homogeneity  and  isotropy  of  the  obtained 
phase  screens  are  quite  good. 

The  behavior  of  the  structure  function  for  different  numbers  of  Zernike  polynomials  used  for  the 
expansion  is  shown  in  Fig.  3.4.  It  demonstrates  that  within  the  confidence  limits  the  structure  functions 
coincide  with  exact  solution.  Note,  that  the  first  two  Zernike  polynomials  make  a  main  contribution  to  the 
phase  fluctuations.  The  contribution  of  higher  order  polynomials,  responsible  for  the  small-scale  fluctuations, 
becomes  significant  at  the  distances  which  are  small  as  compared  to  the  aperture  size.  Fig.  3.5  demonstrates 
that  the  structure  function  approches  exact  solution  as  the  number  of  the  Zernike  polynomials  increases.  Thus, 
the  implementation  of  the  modal  method  with  sufficiently  large  number  of  the  Zernike  polynomials  (=100)  is 
expected  to  allow  consideration  of  the  effects  of  large-scale  as  well  as  small-scale  turbulence  on  statistical 
characteristics  of  the  laser  beam,  propagating  through  the  atmosphere. 

Since  the  structure  function  of  phase  fluctuations  for  the  Kolmogorrov  spectrum  has  a  power 

dependence  on  distance  of  the  type  D—ar^^  where  5=0,10  and  /7=1.67  (for  the  parameters  under 
consideration),  it  is  convenient  to  use  the  logarithmic  scale,  with  the  help  of  which  parameters  a  and  b  may  be 
obtained  from  the  least  square  method.  Fig.  3.6  shows  the  structure  functions  for  different  numbers  of  the 
Zernike  polynomials  at  the  logarithmic  scale.  It  is  seen  that  at  large  distances  (of  the  order  of  grid  size) 
numerical  solution  is  in  good  agreement  with  analytical  for  any  number  of  the  Zernike  polynomials,  while  at 
small  distances  the  satisfactory  agreement  may  be  obtained  only  by  increasing  the  number  of  polynomials.  This 
fact  is  also  confirmed  by  estimating  the  parameters  a  and  b  for  different  number  of  polynomials.  For  two 
Zernike  polynomials  Z?=2.0  and  5=0. 04±  1,  while  for  one  hundred  and  four  Zernike  polynomials  /3=1.78  ±  0,02 
and  ^0.04±0.01.  However,  it  should  be  noticed  that  even  large  number  of  polynomials  does  not  allow  to 
obtain  exact  power  dependence  for  the  structure  function. 

Similar  results  may  be  obtained  for  the  von  Karman  model  of  atmospheric  turbulence  (3.5),  the 

structure  function  for  which  is  defined  by  Eq.(3.16).  Let  the  outer  scale  of  turbulence  be  A=2  m,  the  grid  size 
4,=  Ay  =  A=5A4  cm,  Q  =  10^^^  2=0.53  /mi,  A2=1.2  km,  Ny  =  12S.  The  structure  functions 

obtained  from  Monte  Carlo  simulations  (A#=300  phase  screens)  for  these  parameters  are  shown  in  Fig,  3.7.  For 
comparison  the  Fig.  3.7  shows  the  structure  function  of  phase  fluctuations  on  the  screen  generated  using  the 
FFT-based  method.  It  is  clearly  demonstrated  that  the  Zernike  polynomial  expansion  with  the  large  number  of 
polynomials  gives  better  reproduction  of  the  von  Karman  phase  fluctuations  than  the  FFT-based  method. 

It  may  be  more  emphatically  illustrated  at  the  logarithmic  scale  (Fig.3.8),  that  numerical  solutions  are 
in  good  agreement  with  analytical  solutions  for  different  number  of  polynomials  if  the  distance  ris  of  the  order 
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Fig.  3.3.  Sections  of  the  structure  function  of  phase  fluctuations  obtained  for  J  =  104  Zernike  polynomials  and 
M  =  300  phase  screens.  Curves:  1-  analytical  expression  for  Kolmogorov  model,  2  -  statistical  processing:  2  - 

along  the  x-axis,  3  -  along  the  y-axis,  4  -  along  the  chord  ^  =  /^/3.  Cl  =  5- =  I.O6///72,  A2= 
=7.4km. 


Fig.  3.4.  Examples  of  the  structure  functions  for  different  J.  Curves:  1-  analytical  expression  for  Kolmogorov 
model,  2-5  -  statistical  processing  (M  =  300):  2  -  J  =  5,  3  -  J  =  14,  4  -  J  =  44,  5  -  J  =  104.  Atmospheric 
conditions  are  the  same  as  in  Fig.  3.3. 


Fig.  3.5.  The  structure  functions  at  distances  r  «  A  for  different  J.  Curves  designation  and  atmospheric 
conditions  are  the  same  as  in  Fig.  3.4. 

Fig.  3.6,  The  structure  functions  at  the  logarifmic  scales  for  different  J.  Curves  designation  and  atmospheric 
conditions  are  the  same  as  in  Fig.  3.4. 


Fig.  3.7.  Examples  of  the  structure  functions  for  von  Karman  mode!  of  turbulence.  Curves:  1-  analytical 
expression  for  von  Karman  model,  2-5  -  statistical  processing  (M  =  300):  2  -  J  =  2,  3  -  J  =  14,  4  -  J  =  44,  5  - 

FFT  based  method. =5.44c/72,4  =  2/??,  C]  =  1 A  =  0.53/2/7?, Az  =  \.2knu 

Fig.  3.8.  The  structure  functions  for  von  Karman  model  at  the  logarifmic  scales.  Curves  designation  and 
atmospheric  conditions  are  the  same  as  in  Fig.  3.7. 
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of  aperture  size  A.  However,  for  small  distances  discrepances  may  be  observed  which  indicate  that  small-scale 
fluctuations  are  not  reproduced  quite  adequately. 

In  spite  of  all  difficulties,  the  general  conclusion  can  be  made  that  the  phase  screens  obtained  from 
modal  method  are  applicable  for  the  study  of  laser  beam  propagation  through  the  turbulent  atmosphere.  At  the 
same  time  one  should  remember  that  the  effect  of  the  small-scale  fluctuations  on  the  beam  propagation  will  not 
be  adequately  reproduced. 

4.  Computer  simulations  of  laser  beam  propagation 

In  this  section  we  will  discuss  methodical  aspects  of  numerical  integration  of  the  equation  (2.1) 
describing  the  propagation  of  coherent  light  beam  through  randomly-inhomogeneous  medium.  First,  we  will 
concentrate  on  spectral  method  for  solving  diffraction  equation  (2.12),  (2.14)  in  the  space  between  two  phase 
screens.  Then  we  will  discuss  criteria  of  applicability  of  the  split-step  method.  Finally,  the  method  is 
considered  for  constructing  a  flexible  grid,  which  allows  to  follow  random  wandering  of  the  beam  under  the 
conditions  of  strong  turbulence. 

4.J.  Spectral  meUwd 

Spectral  metliod  based  on  the  fast  Fourier  transform  algorithm  is  currently  one  of  the  most  widespread 
approaches  to  the  problenis  of  quasioptics.  This  comparatively  simple  and  effective  method  has  been 
intensively  discussed  in  the  literature^®-'^-'^.  Nonetheless,  there  exists  a  set  of  methodical  questions,  which  so  far 
have  been  left  aside.  The  objective  of  this  subsection  is  to  clarify  some  of  them. 

Let  us  seek  for  the  solution  to  the  problem  (2.12),  (2.14)  on  the  square  aperture  of  the  size  We 
expand  a  complex  amplitude  of  the  light  field  E(x,y,z)  in  Fourier  series  over  the  square  -4/2£a-^4/2, 
-4i/2^  7^4/2: 


^{px+qy) 

A) 


(4.1) 


Ej{x,y,z)=  2  Z 

p^-Nl2-^\g=-Nn+\ 

When  writing  expansion  (4.1)  we  assume  that  the  function  E(x,y,z)  has  been  replaced  by  a  periodic  function 
with  a  period  4).  Suppose  that  Fourier  image  of  the  function  E  is  either  limited  or  rapidly  decreases  so  that  the 
sum  of  the  first  iV'harmonics  reproduces  E(x,y,z)  wiilx  necessary  accuracy. 

Substituting  (4.1)  into  (2.12)  or  (2.14)  we  have: 


Z2 

2ik  + 

r2n] 

2 

[p"  +9^)4 

exp 

—{px+qy) 

P  Q 

dz 

^  A ) 

A 

From  this  it  follows  that  separate  harmonics  E^^  must  satisfy  the  equations: 


dz 


1 

Ik 


— 

A 


[p^  +q^ 


)K, 


-M2  +  l^ p^NI2,  -N/2  +  l^g£N/2. 
The  solution  to  the  system  (4.3)  has  the  form: 


(4.2) 


(4.3) 


Ki^s.y2)  =  ^  ^  {P'  +9')^^ 


271 


(4.4) 


~Ar/2  +  l^ p^NIl,  --M2  +  l^^^  A^/2. 

ir  the  first  and  the  last  slabs  (i.e.  for  5^=0  and  s=S),  the  step  Az  should  be  replaced  by  the  step  Az/2. 

We  now  introduce  a  two-dimensional  grid  =  wAx,  =  nAx,  where  Ax=  A  I  N'y  0  ^  N-l, 

<n<  N-l.  Then  the  values  of  Fourier  harmonics  E^^(z^_y^^  may  be  computed  from  the  values  of  the 

mplex  amplitude  eJ^x^^^,  y^,  •  Numerically  this  can  be  done  using  the  FFT: 
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J\  a~Q  12=0 

Thus,  knowing  the  field  amplitude  at  z=  ^1/2,  we  may  find  initial  values  of  Fourier  harmonics.  The  next  step 
is  to  compute  the  values  of  harmonics  at  z=  from  equation  (4.4).  The  final  step  is  to  find  the  complex 
amplitude  at  a  distance  z=  zi^m  using  the  inverse  Fourier  transform: 


2ni 

N 


(pm  +  gn) 


MS') 


un 


p^-Nn-^\q^-NI2^\ 
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N 


{pm+qn) 


(4.6) 


In  practice  one  should  remember  that  the  most  standard  subroutines,  computing  the  FFT,  use  the  harmonic 
ordering  from  0  to  N-l.  Since  harmonics  of  the  complex  function  are  periodic,  i.e. 


^N-p,q  “  ^-p,q  ^  ^p,H~q  ^  ^ p-Q  ^  (4.7) 

we  may  easily  rewrite  equations  (4.4)  so  that  they  are  applicable  to  the  calculations  with  the  new  harmonic 
ordering  from  0  to  A^l: 


^p,q  (-^5+1/2)  “  ^p,q  (^5-1/2) 
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for  Nll<  p<  N-l,  0^q<  Nil', 
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4.,(^..i/2)  =  ^^.,(v./2)exp 
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for  NI2<  p<N-l,  NI2<q^N-l. 


The  advantages  of  the  presented  algorithm  are  high  speed  and  the  lack  of  such  typical  for  finite-difference 
schemes  side  effects  as  grid  dispersion  and  diffusion.  Numerical  solution  to  the  diffraction  equation  is  exact  for 

practical  purposes  if  the  function  B  is  infinitely  small  at  the  boundaries  of  the  square  aperture  |  a|  ^  ^  /  2 , 

I  j|  ^  4  /  2  and  its  spectrum  B^  is  infinitely  small  for  |/;|  =  M  2 ,  |  ^  =  M  2 .  When  the  first  condition  is 
violated,  interference  between  initial  beam  and  its  periodic  continuation  occurs.  This  obviously  leads  to  the 
distortion  of  the  beam  profile.  The  violation  of  the  second  condition  causes  spatial  frequency  aliasing,  which  in 

turn  leads  to  the  errors  in  computations  of  the  derivatives  (^B/  and  ^Bl  dy^  in  spectral  domain.  Thus,  in 
order  to  control  the  accuracy  during  the  calculations,  one  should  regularly  estimate  the  values  of  the  function 

1^^  at  the  edges  of  the  grid  as  well  as  the  values  of  its  spectrum  \B^^f  at  Nyquist  frequency. 

4.2.  Criteria  of  the  split-step  method  applicability 

Numerical  analysis  of  the  light  field  propagation  in  the  atmosphere  is  based  on  the  discretization  of 
independent  variables  x,yz  and  representation  of  the  amplitude  and  phase  of  the  field  by  their  values  at  the 
discrete  set  of  points.  Over  the  past  decade  the  question,  whether  this  discrete  model  adequately  describes  the 
propagation  of  a  light  field  in  randomly-inhomogeneous  medium  or  not,  has  been  discussed  in  a  lot  of  papers, 
a  review  of  which  has  been  done  by  Knepp^l  He  formulates  the  following  problems,  which  have  to  be  taken 
into  account  when  choosing  the  number  of  grid  points  and  step  size:  -  accuracy  of  phase  reproduction  on  the 
screen;  -  accuracy  of  the  solution  to  the  wave  propagation  equation;  -  minimization  of  boundaiy  effects. 
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Let  us  discuss  these  questions  in  more  details,  adding  our  own  experience  to  the  ideas,  developed  in 

Ref.  43. 


4.2.1.  Phase  reproduction 

The  choice  of  the  phase  screen  size  and  grid  step  for  the  adequate  reproduction  of  phase  fluctuations 
ranging  between  the  inner  and  the  outer  scales  of  turbulence,  has  been  thoroughly  discussed  in  section  3.  Here 
we  only  note  one  important  condition:  phase  difference  between  the  neighboring  grid  points  must  not  exceed 
the  value  of  rr,  i.e. 


jn,n 

For  the  grid  step  Ax  this  condition  may  be  rewritten  in  the  form: 


AA'max 


dcp 

dx 


<71, 


A;rmax 

X,y 


(4.9) 


(4.10) 


When  generating  the  phase  screen  using  the  spectral  method,  random  field  ^{x,})  is  represented  as  a 
superposition  of  Fourier  harmonics  with  random  amplituides,  the  root-mean-square  values  of  which  are. 


1  I - 2n 

where  ^^(k)  is  the  spectrum  of  phase  fluctuations  on  the  screen,  K^g  =  +  /c/  .  Consequently,  for  the 

harmonic  with  indices  pq 


max 

X 


and  the  condition  (4.10)  takes  the  form: 
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AN 


(4.11) 


A.r<^^[max()c^J/;(K^)r'.  (4.12) 

1  JXJ 

For  the  definiteness,  we  will  use  the  modified  von  Karman  model  to  make  some  estimates.  For  this  model  the 

maximum  of  the  expression  k^F^{k)  occures  at  xtq  =  2;^/  4,  where  4  is  the  outer  scale  of  turbulence.  Taking 
into  account  Eqs.(2.8)  and  (3.6)  we  finally  obtain: 


Ax<3.1 - 4==  (4-13) 

From  this  it  is  clear  that  the  estimate  for  the  grid  step  Ax  depends  not  only  on  the  grid  size  A  and  propagation 
parameters  Jcq,  Q,  lo  but  also  on  the  step  Az in  the  propagation  direction.  This  is  associated  with  the  fact  that 
the  variance  of  phase  fluctuations  on  the  screen  is  defined  by  the  length  of  the  modelled  turbulent  slab. 
Therefore,  additional  ways  of  estimating  Az  should  be  found. 

4.2.2.  Light  wave  propagation 

Since  in  the  space  between  the  screens  light  wave  experiences  natural  diffraction,  the  discretization  of 
z-variable  most  likely  will  not  introduce  essential  error  into  the  simulations,  if  Az  does  not  exceed  the 
diffraction  length  for  the  smallest  inhomogeneity,  reproduced  on  the  screen.  Since  the  size  of  this 
inhomogeneity  is  of  the  order  of  the  grid  step,  a  reliable  estimate  for  Az  has  the  form: 

Az<Jc,Ax\  (4.14) 

Note,  that  this  limitation  is  significantly  stronger  than  that  formulated  in  Ref  43: 

K 

Az  <—Aq Ax 

n  .  (4.15) 

Substituting  estimate  for  Az  from  (4. 14)  into  (4. 13),  we  have 
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For  typical  parameters  of  simulations  (4~1  in,  4  ==0*5  m,  4=5.93-10^  m’^)  we  obtain  that  Ax  must  not 
exceed  several  millimeteres  in  the  case  of  strong  turbulence  ( c/  =  10”^^  )  and  several  centimeters  in  the 

casse  of  weak  turbulence  (C/  =  10'^”^  The  number  of  grid  points  N  must  be  from  several  hundreds  to 

several  tens. 

Another  estimate  for  Ax  may  be  obtained  when  simulating  the  propagation  of  the  focused  beams,  the 
wave  front  of  which  at  z=0  is  given  by: 


q>(x,y)  =  Jc^ 


2  .  2 
X  +y 

2R, 


(4.17) 


where  Rf  is  a  focal  length  of  a  lens.  Taking  into  account  the  condition  (4.9),  we  have: 
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(4.18) 


Ax< - -X- — •. 

■^0  “A)  ^ 

This  condition  may  be  optimized  if  we  replace  by  the  effective  size  of  the  region  of  localization  of  the  beam 
on  the  transmitter.  For  Gaussian  beam  this  size  is  of  the  order  of  10  ao . 


4.2.3.  Boundary  effects 

Numerical  integration  of  the  equation  (2.1),  which  describes  the  propagation  in  free  space,  is 
performed  on  the  grid  of  a  finite  size  4-  To  make  sure  that  numerical  solution  is  valid,  one  should  control  the 
beam  energy  conservation.  This  becomes  crucial  when  the  beam  travels  large  distances  comparable  with  a 

diffraction  length  =  Jcao^ .  The  most  reliable  estimate  of  the  numerical  solution  validity  is  the  analysis  of  the 
beam  spatial  spectrum  while  the  beam  propagates  through  the  phase  screens.  However,  priori  estimates  of  the 
Az  step  are  also  quite  useful. 

Since  the  angle  S  of  the  wave  scattering  and  phase  ^  are  related  by  the  equation 


7  1  5cp  1  5cp 

,  &  = - 0  = - 


(4.19) 


the  energy  of  the  wave,  which  travels  a  distance  Az  in  the  direction  of  propagation,  travels  in  the  perpendicular 
direction  a  distance  SAz.  In  order  to  neglect  the  boundary  effects,  one  should  choose  the  grid  size  not  less 
than  2i9n,axAz.  Using  the  estimate  (4.11)  for  and  assuming  the  modified  von  Karman  model  of  turbulence 
we  have: 


1.334^  4^'^  (4.20) 

It  should  be  noted,  that  the  above  estimate  is  valid  only  for  the  model  of  turbulence  with  limited  outer  scale. 
When  generating  the  phase  screens  using  the  spectral  method,  the  actual  outer  scale  reproducable  on  the  grid 
does  not  exceed  ^.Beam  wander  and  blur  are  not  significant.  In  this  case  boundary  effects  are  negligibly  small 
even  when  Az  is  comparatively  large.  This  fact  is  supported  by  the  Eq.(4.20),  when  substituting  into  it  typical 
parameters  of  simulations.  The  modal  approach  allows  to  use  the  models  with  unlimited  outer  scale  (e.g.  the 
Kolmogorov  and  Tatarskii  models).  In  this  case  estimate  (4.20)  is  not  valid  and  analytical  dependences  for  the 
root-mean-square  shift  of  the  beam  and  effective  radius  should  be  used  to  obtain  the  grid  size  4)  and  step  Az 
(see  section  5). 


4.2.4.  Practical  estimates  of  the  split-step  method  convergence 

In  this  subsection  we  perform  the  numerical  analysis  of  the  split-step  method  convergence.  For 

simplicity,  we  consider  a  one-dimensional  beam  with  complex  amplitude  ^=£(x,z),  the  propagation 
equation  for  which  has  the  form: 


BE  d^E  , 

2*0— =  +  (4-21) 

oz  ax 
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The  spectrum  of  phase  fluctuations  on  a  one-dimensional  phase  screen  may  be  obtained  by  integrating  (3.6) 
over  Ky  and  has  the  form: 

i;(K)  =  exp(-K'  /k^)  (4.22) 

where  Kq  and  k:„  have  been  defined  previously.  For  the  simulations  we  choose  the  following  parameters: 

a^=5ciii,  2.  =  1.06^/2;  4,=409.6cn;  £^=102Acm  I^=5nim  JV=2048,  Ax=2ni, 

Beam  diffraction  length  is  .^=14.82  km,  diffraction  length  for  the  smallest  scale  of  amplitude  fluctuations  is 
a:oAa^=23.7  m. 

From  Monte  Carlo  simulations  we  obtain  the  following  statistical  quantities  describing  the  beam 
propagation: 

the  root-mean-square  displacement  of  the  beam  centroid 

x^={xlY^,  x^  =ixl(x,z)dx/ilix,z)dx,  (4.23) 

long  exposure  effective  beam  radius 

acff  =  (ix‘  <  I{x,  2)>  dxl\ I{x,  z) dx}''^ ;  (4.24) 

coherence  function 

y(x)  =  |(F'(-Ar /  2)£'  (x  1 2))|  /  [<  I{-x  /  2)  ><  I{x  /  2)  ;  (4.25) 

and  the  radius  of  coherence  Pe  defined  at  level  of  y(^)  ; 
the  irradiance  covariance  fluctuations 

b,{x)  =  {{l(,-xl2)-<I{-xl2)>){l{xl2)-<I{xl2)>))l  [O;  {-x  I2)a]{^xl  2)]‘'" ,  (4.26) 

where  aj  =<  /  >  -  <  /  >^  -  is  the  irradiance  variance; 

radius  of  irradiance  covariance  pi  defined  at  level  of  bj{x)\ 
the  log-amplitude  variance  of  at  the  beam  center 

(0)  =<  log*  (.E-(O)  /  E,  (0))  >  -  <  log(jF(0)  /  E,  (0))  >* ,  (4.27) 

where  =  eJ^x,  ^  is  the  field  of  the  unperturbed  beam  propagating  in  vacuum; 
the  scintillation  index  (normalized  irradiance  variance)  at  the  beam  center 

=</(0)>/<7(0)>'  -1,  (4.28) 

All  statistical  quantities  have  been  defined  by  averaging  over  M=100  realizations. 

Figures  4. 1-4.3  show  the  dependences  of  the  defined  quantities  on  the  number  of  phase  screens  S  in  the  strong 
fluctuation  regime:  the  length  of  the  path  z=5  km  (=0.3374  z^),  the  Rytov  variance  at  the  end  of  the  path 
f?Q=\2.11,  From  the  figures  it  is  clear  that  the  energy  characteristics  Xc  and  a^r  have  a  week  dependence  on 
the  number  of  screens  and  converge  rapidly  with  increasing  S.  Note  a  rapid  growth  of  the  light  field 
fluctuations  at  S=32,  when  the  distance  between  the  screens  is  comparable  with  the  diffraction  length  for  the 
inner  scale  of  turbulence.  A  reliable  convergence  of  the  results  reveals  for  S>256,  i.e.  in  the  case,  when  Az 

becomes  less  than  diffraction  length,  defined  from  the  grid  step  (^^  <  ).  The  simulation  results,  given  in 

the  weak  fluctuation  regime  (  z=625  m  (=0.042  confirm  this  conclusion  (see  Figs.  4.4-4.6).  For 

this  case  the  condition  Az=  kl^  yields  S=4,  while  the  condition  Az=  khjd-  yields  S=32.  Thus  the  most  reliable 
estimate  for  the  distance  between  the  screens  is  inequality  (4.14). 

4,3  Flexible  grid 

In  the  strong  fluctuation  regime,  when  the  Kolmogorov  spectral  density  is  chosen,  random  wander  of 
the  beam  in  separate  realizations  may  several  times  exceed  its  initial  radius.  So  in  the  simulations  we  have  to 
use  buffer  zones  in  order  to  minimize  boundary  effects  on  the  grid.  This  results  in  the  increasing  grid  step  and 
decreasing  accuracy  of  computations,  if  the  array  size  remains  unchanged. 

In  this  subsection  we  suggest  a  simple  trick,  that  allows  to  retain  the  high  resolution  of  the  grid  in  the 
strong  fluctuation  regime  without  the  increasing  the  number  of  grid  points.  The  trick  is  based  on  the 
introduction  of  flexible  grid,  which  moves  with  the  wandering  beam.  Below  we  will  discuss  two  ways  of 
handling  the  problem. 
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The  first  way  is  lo  predict  the  movement  of  the  beam  on  the  basis  of  its  wave  front  analysis.  Let  after 
the  passage  of  the  Sth  phase  screen,  the  wavefront  of  the  beam  be  Applying  the  least  square  method  it 

is  easy  to  find  the  tilt  angles  and  along  the  axis  x  a  y  respectively: 


1  £ 


h  zn,tr=l  /t7=l  JCf.  n},o=\  -in=l 


(4.29) 


If  the  beam  propagated  at  angles  und  to  the  initial  optical  axis,  its  center  in  the  (S+l)th  screen 
would  displace  by  a  distance 


(4.30) 


from  its  position  in  the  Sth  screen.  If  we  transform  the  wave  front  uj^x, by  removing  the  tilts 

Us  {x,y)  =  Us  (x,7)  y)  sy^ 

the  beam  will  propagate  along  its  initial  optical  axis.  Repeating  this  procedure  at  each  step  and  remembering 
the  calculated  values  of  the  displacements  (7,)^,  we  obtain  practically  straight-line  propagation  of  the 


beam.  The  algebraic  sum  of  the  values  and  [y,)^  will  give  us  the  displacement  of  the  beam  center  at  the 

end  of  the  path.  The  grid  size  may  be  significantly  decreased  by  reducing  the  buffer  zones. 

The  second  way  is  to  decrease  the  effect  of  accumulation  of  the  displacement,  when  the  beam 

propagates  in  the  medium.  At  each  phase  screen  we  calculate  the  coordinates  of  the  beam  centroid  and 


{7,)  .  Then  the  beam  center  is  shifted  to  the  origin  of  the  coordinates  by  cyclic  reassigning  of  elements  of 

complex  array  As  in  the  first  case,  the  shift  of  the  center  at  the  end  of  the  propagation  path  can  be 

obtained  from  the  sum  of  the  displacements  ^ ,  (7J 

In  our  practice  the  application  of  both  methods  in  the  strong  fluctuation  regime  showed  high 
efficiency. 


5.  The  effect  of  large-scale  refractive  index  fluctuations  on  the  propagation  of 
the  beam  through  the  turbulent  atmosphere 

As  it  was  mentioned  before,  atmospheric  turbulence  has  a  broad  spatial  spectrum  of  the  refractive 
index  fluctuations,  which  cannot  be  completely  reproduced  on  a  grid  of  reasonable  size.  At  the  same  time  this 
spectrum  rapidly  drops  with  increasing  spatial  frequency.  Therefore  the  low-order  aberrations  dominate  in  the 
wave  front  of  the  light  wave.  This  has  been  clearly  demonstrated  in  Ref  45,  where  experimentally  measured 
wave  fronts  of  the  laser  beam,  propagating  in  the  planetary  boundary  layer  are  presented. 

Here  we  employ  the  idea  of  separate  simulations  of  the  large-  and  small-scale  fluctuations  of  the 
refractive  index.  Tliis  approach  enables  us  to  analyze  the  contribution  of  the  different  parts  of  atmospheric 
spectra  on  statistics  of  the  light  beam.  It  is  naturally  to  assume,  that  random  wander  and  blur  of  the  beam  as  a 
whole,  which  are  revealed  in  the  irradiance  profile  during  the  long-exposure  registration,  are  to  a  great  extent 
caused  by  a  low-frequency  part  of  the  spectrum. 

This  section  describes  numerical  simulations  of  these  effects  using  the  modal  representation  of  the 
refractive  index  fluctuations.  The  parameters  of  numerical  algorithm  are  optimized  on  the  basis  of  comparison 
of  Monte  Carlo  simulations  results  with  the  theoretical  predictions.  Various  atmospheric  spectra  and 
fluctuation  regimes  are  considered  and  the  influence  of  the  outer  scales  on  the  energy  characteristics  of  the 
light  beam  is  discussed. 


S.l  Optimization  of  the  phase  screen  model 

In  the  modal  representation  of  atmospheric  inhomogeneities  the  diameter  of  the  circle  D,  over  which 
Zernike  polynomial  expansion  of  a  random  phase  is  performed,  is  a  free  parameter.  At  the  same  time  diameter 


25 


D  defines  the  variances  of  random  Zemike  coefficients,  and  consequently  the  value  of  D  significantly  effects  on 
the  phase  fluctuations  variance  on  the  screen.  It  is  naturally  to  expect  that  in  the  limit,  when  the  number  of 

the  Zemike  modes  tends  to  infinity,  the  dependence  of  cr^  on  D  will  be  decreasing.  Most  likely,  the  ideal  phase 
screen,  wloich  is  reproduced  by  a  sum  of  infinite  number  of  the  Zemike  modes,  has  to  adequately  reproduce  the 
wave  front  at  any  D.  However,  in  simulations  the  number  of  the  Zemike  modes  is  always  limited  and  in  the 
most  cases  not  too  large.  Therefore,  it  is  of  practical  interest  to  estimate  the  diameter  of  the  circle  and  the 
number  of  the  Zemike  polynomials,  for  which  such  effects  like  random  blur  and  wander  of  the  beam,  are 
reproduced  satisfactorily. 

As  a  base  for  the  estimates,  let  us  use  theoretical  predictions,  developed  in  for  the  effective  beam 
parameters  in  the  atmosphere"’*^.  In  particular,  for  the  Kolmogorov  spectral  density,  theory^^  gives  the  following 

formula  for  the  effective  beam  radius  which  defines  the  dimensions  of  the  area  painted  by  beam  wander 

=  +  (5.1) 

where 


r  f 

z 


z 

1- — 


R 


rj 


1/2 


is  the  radius  of  diffraction-limited  Gaussian  beam  in  vacuum,  ao  is  the  initial  beam  radius. 


is  the  Rytov  variance.  In  simulations  we  obtained  the  following  statistical  quantities: 
root-mean-square  displacement  of  the  beam  centroid 

=<  >,  p^  =  =  Wxidxdyl  Widxdy,  7,  =  IJ yldxdy!  Widxdy 

effective  beam  radius 


(5.2) 


(5.3) 


(5.4) 


+7^)  </>  dxdyllhdxd}^^^^ .  (5.5) 

We  performed  the  averaging  over  M=100  realizations  for  the  beam  ^vith  the  following  parameters:  ao=5  cm, 
A=1.06  jj.  m  (diffraction  length  .^=14.82  km).  The  grid  size  was  4)=100  cm  with  128  points  in  both  x-  and  y- 
directions.  The  number  of  phase  screens  was  varied  from  S=10  to  S=20,  the  length  of  the  path  was  z=7.41  km 
(=0.5.^).  We  used  the  Kolmogorov  model  of  turbulence. 


5.1.1.  The  collimated  beam 

When  simulating  the  phase  screens,  we  start  for  simplicity  from  5  Zemike  polynomials,  describeng 
the  abberations  of  the  first  and  the  second  order  (the  piston  mode  is  exclude  from  the  consideration).  Fig.  5.1 

shows  the  dependences  of  on  z,  obtained  for  different  ratios  D!  in  the  case  of  moderate  fluctuation 

regime  (CJ  =5-10"^'^  y^=1.31).  The  best  agreement  between  the  simulation  results  and  theoretical 

prediction  (5. 1)  is  observed  for  D!  ao=2.5.  With  increasing  D  the  effective  radius  decreases  due  to  the  decrease 
in  the  variances  of  the  Zemike  modes,  proportional  to  .  Fig.  5.2  shows  the  same  dependences  obtained  for 
Dl  bq-I.S  in  the  cases  of  weak  (>55=0.26),  moderate  (y55=1.31)  and  strong  (>55=13.1)  fluctuation  regimes. 

It  is  of  interest  to  try  to  improve  the  agreemen  between  the  simulation  results  and  theoretical 
predictions  by  increasing  the  number  of  the  Zemike  polynomials  J.  We  take  J=14  for  the  modes  up  to  the  4th 

order  and  J=44  for  the  modes  up  to  the  6th  order.  The  calculated  dependences  a^{^,  obtained  for  the 
different  number  of  the  Zemike  modes  are  shown  in  Fig.  5.3  (/^=1.31,  Di  a^-A)  and  in  Fig.  5.4  (>C^=13.1, 

D!  ^0=8).  From  these  figures  it  is  clear  that  increasing  the  number  of  the  Zemike  modes  gives  the  more 
pronounced  effect  in  the  case  of  moderate  fluctuations  and  less  pronounced  in  the  case  of  strong  fluctuations. 
Note,  that  the  optimum  size  of  the  diameter  D  increases  with  increasing  number  of  the  Zemike  polynomials 
from  D!  ^0=2.5  for  J=5  to  Di  ao=S  for  J=44. 
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Fig.5.1.  Fig. 5.2. 

Fig. 5.1.  Dependencies  of  on  Z  for  collimated  beam.  Curves:  1  -  theory,  2-4  -  Monte  Carlo  method  at 
different  values  of  D!  2  -  D!  a^^2.5\  3  -  D!  a^^4\  4  -  D!  a^=10\  5  -  vacuum.  Turbulence  is  reproduced  by 
J=5  Zernike  modes.  Propagation  conditions  are:  a^  =Scm,X  =  \.06}4m,C^  =  5-10"*’ =  1.31. 

Fig. 5.2.  Dependencies  of  a^^  on  Z  for  collimated  beam.  Curves:  1,2,3  -  theory,  T,2',3‘  -  Monte  Carlo  method 
((^  =5•10’'^5*10“^^,10“^^c/n~^^^  respectively).  Phase  screen  parameters:  J=5,  Di  5q=2.5.  Beam  parameters: 
aQ-5cm,X  =  I.O6///27. 


Fig.5.3,  Fig.5.4. 

Fig. 5.3.  Dependencies  of  a^^  on  Z  for  collimated  beam.  Curves:  1  -  theory,  2-4  -  Monte  Carlo  method  at 
different  number  of  Zernike  modes:  2  -  J=44:  3  -  J=74;  4  -  J=5;  5  -  vacuum.  Phase  screen  parameter  Di  ^0=4. 
Propagation  conditions  are:  a^  =  ‘Scn\X  =  1.06/i/72,  =5-10^*^  ,0^q  =  1.31. 

Fig.5.4.  Dependencies  of  a^^  on  Z  for  collimated  beam.  Curves  designation  is  the  same  as  in  Fig.  5.3.  Phase 
screen  parameter  D!  a^-2>.  Propagation  conditions  are:  a^  =  5c/72,A  =  1.06/;/j3,  =  5*  10“^^  =  1.31. 


Fig. 5.5.  Fig. 5.6. 

Fig. 5. 5.  Mean  irradiance  profiles  obtained  by  Monte-Carlo  method:  1  -J=5;  2  -  J='/4:  3  -  J=44  (Df  3^=4).  Curve 

4  -  Gaussian  profile  with  theoretically  predicted  a^^. 

Fig. 5.6.  Logarithmic  representation  of  mean  irradiance  profiles  shown  in  Fig. 5. 5. 
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Drastic  discrepances  between  numerically  obtained  and  theoretically  predicted  resuilts  in  the  regime 
of  strong  fluctuations  (Fig.  5.4)  may  not  be  explained  by  the  limited  number  of  the  polynomials:  the  increase  in 

J  from  14  to  44  causes  only  2%  increase  in  From  our  point  of  view,  the  theory  overestimates  the  values  of 


for  the  conditions  under  discussion.  Indeed,  the  results  in  Ref  46  are  obtained  under  the  assumption  that 
the  long-exposure  mean  irradiance  profile  is  Gaussian.  Analysis  of  the  mean  irradiance  profiles  obtained  in 
simulations  for  different  J  (Fig.  5.5)  demonstrates  that  in  the  strong  fluctuation  regime  this  assumption  might 
be  not  valid.  This  may  be  clearly  demonstrated  by  the  analysis  of  the  Fig.5.6,  where  the  functions 


are  shown  for  the  mean  irradiance  profiles  from  Fig.  5.5.  Here  is  the  maximum 


value  of  the  intensity  for  the  given  profile.  A  straight  line  on  the  plot  corresponds  to  Gaussian  profile  with  the 
effective  radius  obtained  from  Eq.(5.1).  The  difference  between  the  slopes  of  the  simulated  and  theoretically 
predicted  lines  is  clearly  seen.  This  diifference  is  the  most  pronounced  at  the  edge  of  the  beam. 

Bearing  in  mind  the  considerations  presented  above,  it  is  likely  to  expect  that  with  increasing 
fluctuation  strength  the  difference  between  the  simulated  and  theoretically  predicted  (5.1)  radius  of  the 
effective  beam  will  be  increasing. 


5.1.2.  Focused  beam 

When  analyzing  the  propagation  of  the  focused  beam,  we  use  the  simplified  phase  screen  model, 
which  contains  only  the  first  and  the  second-order  Zernike  modes  (J=5).  The  optimized  ratio  D!  5o  for  this 
case  is  equal  to  2.5.  Let  us  consider  the  beam,  described  in  subsubsection  5.1.1,  focused  at  a  distance  Rf-lAl 

km  (=0.5.^).  Simulated  and  theoretical  dependences  of  on  z  are  shown  in  Fig.  5.7  for  those  values  of  6^, 

which  correspond  to  /55=0.26,  1.31,  13.1.  It  is  seen  that  for  the  focused  beam  as  well  as  for  the  collimated 
beam,  numerical  and  theoretical  results  quite  satisfactorily  agree  in  the  weak  and  moderate  fluctuation 

regimes,  while  in  the  strong  fluctuation  regime  the  formula  (5.1)  overestimates  a^.  At  the  same  time  the 
location  of  the  effective  beam  waist  in  turbulence  obtained  from  the  Monte  Carlo  simulations  is  in  good 
agreement  with  that  predicted  by  Eq.(5.1)  in  wide  range  of  CJ  (Fig.5.8).  As  predicted  in  Ref.46,  turbulent 
beam  waist  migrates  towards  the  transmitter  as  the  turbulent  strength  increases. 

5.2.  The  role  of  the  outer  scale 

The  estimate  of  the  influence  of  the  outer  scale  of  turbulence  on  energy  characteristics  of  the  beam 
was  performed  for  the  von  Karman  spectral  density.  (3.6).  The  conditions  of  the  beam  propagation  are  those 
described  in  subsection  5.1. 

Fig.  5.9  shows  the  simulated  dependences  obtained  for  the  collimated  beam  in  the  moderate 

fluctuation  regime  (Q  =5-10"*’  cm’^/^)  for  different  values  of  the  outer  scale  4.  For  comparison,  we  again 

plotted  the  curves  for  the  Kolmogorov  turbulence  presented  in  Fig.  5.1:  simulated  dependence  for  4  = 
theoretical  prediction  (5.1).  It  is  clear  that  with  decreasing  outer  scale  the  effective  beam  radius  decreases.  This 

decrease  is  more  pronounced  for  small  4  •  Basically,  this  behavior  of  a^^  can  be  explained  by  the  decrease  in 
the  variance  of  random  migrations  of  the  beam  and  to  a  less  degree  by  the  decrease  in  random  blur,  since  the 
value  of  the  outer  scale  most  significantly  affects  the  tilt  variance  of  the  wave  front^^*'’^ 

Results  for  the  focused  beam  are  presented  in  Figs. 5. 10,  5.11,  where  7?/=0.5.^.  Fig.  5.10  shows  the 
dependence  in  the  weak  (>^=0.26)  and  moderate  (>62=1.3 1)  fluctuation  regimes  for  both  the 

Kolmogorov  and  the  von  Karman  (with  4  =  100^0=5  m)  models  of  turbulence.  For  comparison  the  curves 
obtained  from  Eq.(5.1)  are  also  presented.  Fig.  5.11  illustrates  the  dependence  of  the  location  of  the  effective 
beam  waist  zic  on  the  outer  scale.  Results  were  obtained  from  Monte  Carlo  simulations  in  the  regimes  of  weak, 
moderate  and  strong  fluctuations.  With  decreasing  outer  scale,  the  effective  beam  waist  shifts  from  the 
transmitter  in  all  regimes  according  to  the  Kolmogorov  model. 
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Fig.5.7.  F\g.5.Q. 

Fig. 5.7. Dependencies  of  on  z  for  convergent  beam.  Curves:  1,2,3  -  theory,  1',2',3'  -  Monte  Carlo  method 

(/^  =  13.1,1.31,0.26,  respectively).  Phase  screen  parameters:  J=5,  Dl  aQ=2.5.  Beam  parameters: 
Sq  =5cni,X  =  1.06/7/72,  Rj.-lA\kin(=0.5zj), 

Fig. 5. 8.  Coordinate  of  effective  beam  waist  in  Kolmogorov  turbulence  a  s  a  function  of  6^.  Curves:  1  -  theory,  2 
-  Monte  Carlo  method.  Phase  screen  parameters:  J=5,  Dl  aQ=2.6.  Beam  parameters: 
/7q  =  Scm,  X  =  1.06/7/72,  Rjr  =  7.41.to(=  0.5z^). 


Fig.5.9.  Fig. 5.10. 

Fig. 5.9.  Dependencies  of  on  Z  for  collimated  beam.  Curves:  1  ,2-  theory  and  Monte  Carlo  method  for 
To  =00.  3,4  -  Monte  Carlo  method  for  ^  =5m  and  2  m  respectively,  5  -  vacuum.  Propagation  conditions  are: 
=  5c/7^A  =  1.06/7/7;  =  5*10“^’ c/77-''\/^  =  1.31. 

Fig. 5.10.  Dependencies  of  on  Z  for  collimated  beam.  Curves:  1  ,2-  theory  and  Monte  Carlo  method  for 
=00,  =  5*  10''^^c/72"^^^,  3  -  Monte  Carlo  method  for  =5/77,  Cj  =  5- 10“^^ c/72~^^^,  4,5  -  theory  and  Monte 

Carlo  method  for  ^  =oo,  Cj  =  1 0"*^  c/77■^^^  6  -  Monte  Carlo  method  for  =  5/72,  =  10■^^c/72“^^^  7  -  vacuum. 

Propagation  conditions  are:  a^  -  5cm,  A  =  I.O6/7/72,  R^.  =  7.4\km{-  0.5z^) . 


Fig. 5.11. 


Fig.5.11.  Coordinate  of  effective  beam  waist  in  von  Karman  turbulence  a  s  a  function  of  Curves:  1  - 
=  10■^^c/72■^^^4  =  0*26.  2  -  =5•10■^’c/72■^^^^  =  1.31,  3-  6?  =  5-10"^^  c/77"^^\y^  =  13.1.  Dashed  lines- 

values  of  I^=  00,  Beam  parameters:  a^  =  Scin,X  =  I.O6/7/72,  R^^lA\kw{=  0‘5z^) 
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S.3.  Random  migrations  of  beam  centroid 

Monte  Carlo  method  allows  to  separately  investigate  various  statistical  quantities  describing  the  beam 
propagation.  In  particular,  this  method  may  be  used  to  confirm  the  results,  predicted  in  Ref.  50  for  the  variance 
of  random  migrations  of  the  beam  in  the  turbulent  atmosphere.  Under  the  assumtion  of  the  disturbed  field  of 
Gaussian  collimated  beam  authors  of  Ref.50  derived  the  formula  for  the  variance  of  the  beam  centroid 
displacement: 


^23 

471  z 


(k)^  ck. 


(5,6) 


where  is  the  spectrum  of  the  refractive  index  fluctuations.  For  the  Kolmogorov  model  of  turbulence  the 
integral  in  (5.6)  may  be  computed  analytically.  Thus  we  have 


,  0.132u^r(l/6)  ^  3  .,,3 

Pc  3  2^'^  ' 


(5.7) 


For  other  than  Kolmogorov  spectral  densities,  integral  in  Eq.(5.6)  may  be  easily  computed  numerically. 

The  Table  5.1  presents  the  results  of  numerical  analysis  of  migration  of  collimated  beam  centroid 
random.  The  data  are  compared  with  those  obtained  from  Eqs.(5.6),  (5.7).  Good  agreement  between  the 
theoretical  predictions  and  results  of  Monte  Carlio  simulations  is  demonstrated  in  the  wide  range  of  values  of 
the  outer  scale  and  turbulence  strength. 


Tabie  5.1,  RMS  beam  displacement  as  a  function  of  outer  scale  of  turbulence  for  collimated  beam 
(5o=5c/2^  /I  =  1.06/7/2^  z=l  Alknf^O.Sz^)), 


4  /  3(1  =  CO 

Theory  Monte  Carlo 

4/a„=100 

Theoiy  Monte  Carlo 

A/ao=40 

Theory  Monte  Carlo 
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0.356 

0.299 

0.314 

0.240 

0.244 

5-10“'' 

1.31 

0.800 

0.806 
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0.708 

0.53 

0.55 

5- 10'“ 

13.1 

2.58 

2.55 

2.11 

2.18 

1.77 

1.75 

6.  Small-scale  turbulence 


Fig.  6.1,  Spectra  of  the  refractive  index 
fluctuations  in  the  atmosphere.  A  -  the  Andrews 
spectrum.  K  -  the  modified  von  Karman 
spectrum. 


In  this  section  we  analyze  the  effect  of  small-scale 
fluctuations  of  the  refractive  index  on  statistical  characteristics  of 
the  laser  beam  propagating  in  the  turbulent  atmosphere.  We 
consider  the  changes  in  the  following  characteristics  with  the 
distance:  the  coherence  of  light  field,  the  profile  of  the  mean 
irradiance,  irradiance  variance  and  covariance  in  case  of  the 
collimated  and  focused  beams  at  different  parameters  of 
turbulence.  The  results  of  Monte  Carlo  simulations  are  compared 
with  predictions  of  the  scintillation  theory^  and  phase 
approximation  in  Huygens-Kirchhoffs  method^. 

The  obtained  statistical  characteristics  are  the  results  of 
processing  of  the  ensemble  of  light  field  distributions  in  the 
registration  plane.  In  calculations  for  each  realization  we  follow 
a  random  blurring  of  the  beam  and  put  its  axis  at  the  origin  of 
the  coordinates. 

To  perform  the  analysis  of  the  effect  of  different  models  of 
atmospheric  turbulence  on  statistical  characteristics  of  the  beam 

we  considered  the  modified  von  Karman  O^'^(k)  (3,6)  and  the 

Andrews  (3.7)  spectra.  For  both  models  the  phase 

screens  were  generated  using  the  same  arrays  of  random 
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numbers  in  order  to  exclude  random  deviations  associated  with  a  finite  number  of  realizations  M 

The  presented  simulations  are  made  for  initially  Gaussian  beams  and  laser  wavelength  of  0.5  /rni. 

Computer  simulations  parameters 

•  In  the  case  of  the  collimated  beam  : 


initial  beam  radius 

^0=2  cm 

length  of  the  path 

2=2  km 

inner  scale  of  turbulence 

/o  =4  mm 

outer  scale  of  turbulence 

4=25  cm 

structure  constant 

the  irradiance  variance 

calculated  from  the  Rytov  theory  for  2=2  km 

/?o=17.1 

structure  function  of  phase  fluctuation  for  spherical  wave 

on  the  transmitting  aperture  for  z=l  km 

A(2a)=104.6 

Grid  parameters: 

number  of  screens 

S=1Q 

distance  between  the  screen 

A2=100  m 

the  variance  of  phase  fluctuations  on  the  screen 

c4=3.83  (0/“'(k)) 

grid  spacing 

ol=4.U 
h=l.O  mm 

number  of  grid  points 

74 X  74=512x512 

size  of  the  grid 

A=  Ak=51.2  cm 

number  of  realisations 

A^80 

For  the  chosen  grid  spacing  the  smallest  scale,  that  can  be  represented  on  the  grid,  is  ro  -2mm.  This 
provides  an  adequate  reproduction  of  high-frequency  structure  of  fluctuations  in  the  modified  von  Karman 

and  Andrews  spectra,  which  are  shown  in  Fig.  6.1  at  the  logarithmic  scale. 

•  In  the  case  of  the  focused  beam : 


initial  beam  radius 

ao=2  cm 

focal  length  of  the  lens 

Z/=500  m 

length  of  the  path 

2=1  km 

inner  scale  of  turbulence 

Iq  =4  mm 

outer  scale  of  turbulence 

4=13  cm 

staicture  constant 

the  irradiance  variance 

Cj,  =  5-10*'^  cm 

calculated  from  the  Rytov  theory  for  2=1  km 
structure  function  of  phase  fluctuation  for  spherical  wave 

y0§=8.O3 

on  the  transmitting  aperture  for  2=1  km 

A(2a)=87.2 

Grid  parameters: 

number  of  screens 

5^20 

distance  between  the  screen 

A2=50  m 

the  variance  of  phase  fluctuations  on  the  screen 

a^^=1.02  (0/*(k)) 
c4=1.15  i<^/{k)) 
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grid  spacing 
number  of  grid  points 

size  of  the  grid 
number  of  realisations 


^0.25  mm 
A^^xiy,=512x512 

4^  =  A=12.8  cm 
71^120 


In  the  case  of  the  focused  beam  the  grid  spacing  of  0.25  mm  allows  to  reproduce  the  beam  profile  in  the 
focal  plane  of  the  lens.  In  vacuum  the  radius  of  the  focal  spot  is  0.2  mm. 

When  the  grid  spacing  Az  along  the  propagation  direction  is  constant,  the  number  of  screens  S 
decreases  as  the  receiver  plane  approaches  towards  the  transmitting  aperture.  However,  for  5>  3^5,  the  phase 
screen  model,  according  to  the  section  4  of  the  present  paper,  is  valid. 


6J.  Field  renlizations 

The  typical  distribution  of  the  random  irradiance  in  the  beam  cross-section  at  dilferent  distances  from 
the  transmitting  aperture  is  shown  in  Figs.  6.2,  6.3.  The  distortions  in  the  beam  profile  such  as,  for  example, 
local  random  focusing,  caused  by  small-scale  fluctuations  of  the  refractive  index,  are  clearly  demonstrated. 

For  the  collimated  beam  the  effect  of  the  "bump"  in  the  Andrews  spectrum  (3.7)  is  observed.  At  the  end 
of  the  propagation  path  the  beam  breaks  up  into  small-scale  speckles  (Fig.  6.2b).  In  the  case  of  the  modified 
von  Karman  spectrum  (3.6)  the  irradiance  distribution  is  smoother  and  the  characteristic  scale  of 
inhomogeneities  is  larger  than  in  the  case  of  the  Andrews  spectrum. 

In  the  case  of  the  focused  beams  the  difference  between  the  irradiance  distributions  obtained  for 

different  atmospheric  spectra  is  not  so  distinct  (Fig.6.3)  since  strong  focusing  (7^//  4^=0. 1)  to  a  great  extent 
suppresses  phase  fluctuations. 


6.2.  Coherence  of  light  Jleld 

Spatial  coherence  of  the  light  field  is  characterized  by  the  function  y(p).  At  the  beam  center  this 
function  is  given  by: 


{E{-pn)E-{pii)) 
({l{-pl2)){l{pl2))Y 
where  (  )  denote  averaging  over  an  ensemble  of  Af  realizations. 


Fig.  6.4.  The  coherence  function  |y(/))|  of  the  beann  at  different  distances  for  the  transmitter. 

a)  collimated  beam 

b)  focused  beam 
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100  m 


500  m 


Fig.  6,2a 

Fig.  6.2a  The  propagation  of  one  realization  of  the  collimated  beam  through  the  modified  von  Karman 
turbulence.  Irradiance  distributions  are  registered  in  different  z-planes  along  the  propagation  direction.  Irradiance 
is  normalized  to  a  maximum  value  1^  in  the  planes: 


z,  m 

LJh 


100  500  1000  1500 

1.3  3.3  1.6  1.0 
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z=  100  m 


z  =  500  m 


Fig.  6.2b 

Fig.  6.2b  The  propagation  of  one  realization  of  the  collimated  beam  through  the  Andrews  turbulence.  Irradiance 

distributions  are  registered  In  different  z-planes  along  the  propagation  direction.  Irradiance  is  normalized  to  a 

maximum  value  I,„  in  the  planes: 

z,  m  100  500  1000  1500 

4//0  1.3  3.6  1.2  0.7 
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Fig.  6.3a 

Fig.  6.3a  The  propagation  of  one  realization  of  the  focused  beam  through  the  modified  von  Karman  turbulence. 
Irradiance  distributions  are  registered  in  different  z-planes  along  the  propagation  direction.  Irradiance  is 
normalized  to  a  maximum  value  /„  in  the  planes: 


Z  =  50  m  z  =  250  m 


Fig.  6.3b 

Fig.  6.3b  The  propagation  of  one  realization  of  the  focused  beam  through  the  Andrews  turbulence.  Irradiance 
distributions  are  registered  In  different  z-planes  along  the  propagation  direction.  Irradiance  is  normalized  to  a 
maximum  value  in  the  planes: 
z,  m  50  250  500  750 

IJh  2.0  12.9  22.0  5.6 
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Fig.  6.4  shows  the  module  of  the  function  j^(p)  in  the  cases  of  collimated  and  focused  beam.  The  curves 
are  obtained  from  polynomial  interpolation  of  experimental  data: 

/(p)=l  +  ^a,p?''  (6-2) 

1=1 

In  the  case  of  the  collimated  beam  the  function  \y{p)\  at  small  distances  3  mm  becomes  narrower  as 

the  light  propagates  along  the  z-axis  (Fig.  6.4a).  At  the  beginning  of  the  path  (z=550  m)  the  function  |y(p)|  is 
significantly  narrower  for  the  Andrews  spectrum  than  for  the  modified  von  Karman  spectrum.  The  reason  for 

this  is  the  larger  contribution  of  small-scale  fluctuations  to  the  spectrum  (3.7). 

For  2=1050  m  the  functions  |y(p)|  virtually  coincide  for  both  models  of  atmospheric  turbulence.  At  the 

end  of  the  propagation  path  (2=1550  m)  the  difference  between  the  functions  |7{p)|  for  the  models  under 
discussion  is  revealed  only  at  distances  p>2  mm.  For  p>4  mm  and  2=1550  m  the  rate  of  the  decrease  in 
coherence  slows  down  due  to  natural  divergence  of  the  beam. 

In  the  focused  beam  the  changes  in  |^(p)|  with  z  are  more  distinct  (Fig.  6.4b).  Up  to  the  focal  plane 

7?/=500  the  width  of  the  function  |y(p)|  decreases  due  to  the  focusing  and  atmospheric  turbulence.  The 
broadening  of  the  coherency  function  after  the  focal  plane  is  caused  by  a  divergence  of  the  beam  as  a  whole. 
The  effect  of  "bump"  in  the  Andrews  spectrum  is  revealed  at  distances  z<Rf.  In  this  case  the  width  of  the 
coherency  function  is  less  than  in  the  case  of  the  von  Karman  spectrum.  When  ^  Rf  the  discrepancies  between 
coherency  functions  obtained  for  two  models  of  turbulence  are  not  large. 

In  the  simulations  |;'(p)|  changes  from  1  down  to  0.6-0.7.  To  estimate  the  radius  of  coherence  pe  at 
small  distances  p ,  we  will  use  a  parabolic  approximation  of  |y(p)|  in  the  vicinity  of  p  =0: 

|r(p)|  =  i+V- 

Defining  the  radius  of  coherence  at  e''  level  of  the  function  \r{j^\ 

7Wl=«-'. 

we  obtain  the  following  expression  for  p^'. 


The  change  of  the  radius  of  coherence  pg  with  z  for  the  collimated  and  focused  beams  is  shown  in 
Fig.(6.5). 

In  the  collimated  beam  the  radius  ps  rapidly  falls  at  a  distance  z<2'=800  m.  Its  further  changes  with  the 

distance  are  insignificant  (Fig.  6.5a).  At  a  distance  z=z‘  the  radius  of  the  first  Fresnel  zone  Zf  =  ^[II  is  2.0 
cm.  This  value  coincides  with  the  initial  beam  radius  ao.  Thus,  at  a  distance  z  the  decrease  in  correlation  of 
the  light  field,  caused  by  joint  effects  of  the  turbulence  and  diffraction,  stops. 
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Fig.  6.5.  Radius  of  coherence  ps  as  a  function  of  z.  Symbol  "K"  denotes  for  the  modified  von  Karman 
turbulence  and  "A"  for  the  Andrews  turbulence. 

a)  collimated  beam; 

b)  focused  beam 


At  the  very  beginning  of  the  propagation  path  zj^150  m  the  radius  Pe  for  the  collimated  beam  may  be 
estimated  from  the  plane  wave  approximation.  Fig.  6.5a  shows  analytical  dependence  of  the  radius  of 
coherence  on  z  for  the  case  of  the  plane  wave  and  the  Kolmogorov  spectrum  of  atmospheric  turbulence^: 


9  EP  ~ 


1 

{lASClk^zf^ 


(6.6) 


For  2<800  m  (z<  0.15Zt/)  the  radius  of  coherence  pe  of  the  beam  is  close  to  that  of  the  plane  wave  pEp.  For 
2>800  m  the  natural  divergence  of  the  beam  compensates  for  the  decrease  in  correlation  in  the  turbulent 
atmosphere,  and  the  radius  Pe  remains  practically  unchanged.  As  a  result,  for  2>800  m  radius  of  coherence  Pe 
of  the  beam  is  larger  than  the  radius  of  coherence  Pep  of  the  plane  wave.  The  lower  boundary  of  the  region  ^ , 
where  ps  >  Pep,  may  be  estimated  from  the  condition^: 


(6.7) 

z 

For  the  parameters  under  discussion  /  »700  m,  that  is  close  to  the  result  of  numerical  simulations.  However, 

this  value  of  z  more  likely  satisfies  the  weak  inequality  —  <  1.22pQ^^^  rather  than  the  strong  inequality  (6.7). 

z 


6,3  Distribution  of  the  mean  irradlance 

Mean  irradiance  profiles  were  obtained  for  a=0.  For  the  collimated  beam  the  irradiance  (/(p))  at 

different  distances  z  is  shown  in  Fig.  6.6.  Although  the  number  of  realizations  is  quite  large  (A#=80), 
irradiance  profile  reveals  several  peaks,  associated  with  the  error  of  Monte  Carlo  simulations. 

The  effect  of  the  small-scale  turbulence  on  the  distribution  of  irradiance  in  the  beam  is  characterized  by 
a  mean  short-exposure  radius  a.  In  simulations  is  computed  from  Eq.(5.1). 

The  radius  a  is  an  integral  parameter,  that  smoothes  the  irradiance  over  the  beam  section.  For  the 
collimated  beam  the  value  a'^  steadily  grows  with  z  (Fig.  6.7a).  Diffractive  broadening  of  the  radius  a  is  not 
large:  it  forms  approximately  7%  for  z==1.9  km.  Due  to  the  turbulence  the  beam  broadens  by  a  factor  of  3-4 
while  it  propagates  along  the  same  length.  The  Andrews  model  provides  the  larger  broadening  than  the 

modified  von  Karman  model.  Fig.  6.7a  shows  analytical  estimate  (5.1)  of  obtained  from  the  theory  of 
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Fig.  6.6.  Mean  irradiance  {l{p))  of  the  collimated  beam  in  the  plane  x=0.  a)  in  the  modified  von  Karman 
turbulence;  b)  in  the  Andrews  turbulence. 
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strong  fluctuations  in  Ref.  46,  The  analytical  estimate  of  that  includes  both  random  motion  blurring  of  the 

beam,  exceeds  the  result  of  simulations  in  the  conditions  of  small-scale  turbulence. 

The  obtained  from  Monte  Carlo  simulations  dependence  of  a  on  distance  z  enable  us  to  find  the  mean 
radius  of  curvature  introduced  in  Ref  46.  Indeed, 


da{z)ldz 


(6.8) 


For  the  collimated  beam  the  change  of  with  distance  is  shown  in  Fig.  6.8a.  When  z  >500  m  the 
dependencies  for  the  Andrews  and  von  Karman  spectrum  coincide  with  theoretical  prediction"^.  A  discrepancy 
between  theoretical  and  simulated  values  of  for  z  <  500  m  arises  since  our  numerical  algorithm  does  not 
lake  into  consideration  random  blurring  of  the  beam,  though  it  is  precisely  this  effect,  that  is  responsible  for 
the  change  in  R,  at  the  beginning  of  the  propagation  path.  The  discrepancy  between  the  values  of  Rc  at  the 
beginning  of  propagation  explain  the  decrease  in  mean  radius  a,  as  compared  with  theoretical  predictions,  at 
the  end  of  the  path  (Fig.  6,7a). 

For  the  focused  beam  simulated  values  of  R  coincide  with  theoretical  valucs"^^  for  z<200  m  z>500  ni 
(Fig.  6.8b).  Some  discrepancies  in  the  vicinity  of  the  beam  waist  are  most  likely  associated  with  the  error  in 
numerical  algorithm  and  require  further  analysis. 


The  coordinate  of  the  effective  beam  waist  may  be  computed  from  the  equation 


=  0 ,  For  both 


dz 

models  of  turbulence  the  beam  waist  shifts  towards  the  transmitter.  The  simulated  shifts  are  less  than 
theoretically  predicted.  The  reason  for  this  is  underestimating  the  role  of  large-scale  fluctuations  in  numerical 
algorithm. 


In  the  case  of  the  focused  beam  the  peaks  of  the  mean  irradiance  (/(>c>))  are  not  high,  since  for 


Rf  i  Zty=0. 1  the  behavior  of  the  beam  is  governed  by  the  focusing.  Note,  that  in  the  focal  plane  peak  value  of 
the  mean  irradiance  is  significantly  larger  for  the  modified  von  Karman  model,  than  for  the  Andrews  model. 
The  dependence  of  a'^  on  distance  for  the  focused  beam  is  plotted  in  Fig.  6,7b.  As  in  the  case  of  the 
collimated  beam,  theoretical  estimate  of  a  in  the  case  of  the  focused  beam  exceeds  the  values  obtained  from 
Monte  Carlo  simulations. 
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Fig.  6.9.  Mean  irradiance  the  focused  beam  in  the  plane  x=0.  a)  in  the  Andrews  turbulence;  b)  in  the 

modified  von  Karman  turbulence. 


6,4.  Irradiance  variance 

The  normalized  irradiance  variance  is  computed  from  the  equation: 
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(6.9) 


y9/(A^) 


{l\p,z))-{l(p,z)y 


^  {I(p,z)f 

For  simplicity  of  notation,  we  will  suppress  the  z  dependence  unless  helpful  for  clarity.  In  simulations  we 
consider  irradiance  variance  at  the  beam  center  p  =  0  and  at  the  diffractive  beam  edge  p- For  the 
collimated  beam  the  normalized  irradiance  variance  (?,  grows  at  the  beginning  of  the  propagation  path  and 
saturates  when  ir  >1000  m  (Fig.  6.10a).  The  normalized  irradiance  variance  f?,  is  always  larger  at  the 
diffractive  beam  edge  p=  than  at  the  beam  center.  This  is  in  complete  agreement  with  experimental  data 
and  theoretical  estimates^.  Note,  that  for  both  the  Andrews  and  the  von  Karman  models  of  turbulence  the 
dependencies  ^(z)  are  nearly  identical. 


Fig.  6.10a  Fig.  6.10b 


Fig.  6.10.  Irra(diance  variance  for  (a)  collimated  and  (b)  focused  beams. 

In  the  case  of  the  focused  beam  the  normalized  irradiance  variance  fi]  grows  only  at  the  very  beginning 
of  the  propagation  path  (z  <Rf)  (Fig.  6.10b),  In  the  vicinity  of  the  effective  beam  waist  in  the  turbulent 

atmosphere  at  the  beam  center  and  at  the  diffractive  beam  edge  some  decrease  in  is  observed. 

Past  the  beam  waist,  the  irradiance  variance  at  the  beam  center  Pl{o)  slowly  increases  with  the 

distance,  while  the  variance  at  the  diJQfractive  beam  edge  P]  {a^)  climbs  sharply. 

Thus,  we  come  to  a  conclusion  that  field  of  irradiance  fluctuations  in  the  section  of  the  focused  beam  is 
close  to  statistically  homogeneous  in  the  region  of  weak  fluctuations  (P^  <1)  and  in  the  small  region  beyond  the 
effective  beam  waist  2i,c.  Immediately  in  the  beam  waist  z«  25,^  and  far  beyond  it  z>  4^,  region  of 

strong  fluctuations  (y?o^»l),  the  field  of  irradiance  fluctuations  in  the  beam  section  is  statistically 
inhomogeneous.  The  latter  contradicts  the  results  obtained  from  the  phase  approximation  in  Huygens- 
Kirchhoffs  method^  and  the  scintillation  theory  based  on  the  effective  beam  parameters'^^,  further  referred  as 
the  "effective"  theory. 

Monte  Carlo  method  and  the  scintillation  theory  of  the  effective  beam  parameters 
To  compare  the  results  of  Monte  Carlo  simulations  with  results  of  the  scintillation  theory  of  the  effective 

beam  parameters,  let  us  consider  the  behavior  of  the  log-irradiance  variance  aj^^(p,z)  given  by: 
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"0  J 


h  j 


(6.10) 


The  log-irradiance  variance  has  been  calculated  at  the  beam  center  p=0  and  at  the  diffractive  beam  edge 


Following  Miller  et.  al"^,  let  us  first  analyze  the  case  of 
the  weak  turbulence  when  =  5- 10"‘^  cm■^^^  Fig.  6.11 
shows  the  log-irradiance  variance  at  the  beam  center 
/?=0  and  at  the  diffractive  beam  edge  p-  Bd  of  the  focused 
beam  as  a  function  of  parameter  O  =  zl  ka\ ,  for  parameter 
klf-  Rflka]  =0.1.  The  variance  begins  by  closely 

following  the  Rytov  variance  Pl  in  the  region  Q<Q/, 
plunging  slightly  to  a  minimum  in  the  near  vicinity  of 

0.  =  klf.  Past  Q.-D.fy  the  log-irradiance  variance 
climbs  sharply  to  a  local  maximum,  then  grows,  all  the  way 
remaining  less  than  the  Rytov  variance  pl .  At  the  diffractive 
beam  edge,  the  climbs  more  rapidly  when  Q<Q/. 

When  Q  is  near  Q/  the  variance  climbs  sharply  to  a 

local  maximum  and  then  drops  at  Q  =  Q/. 

In  simulations  the  end  of  the  propagation  path  2=1  km 
corresponds  to  f2=0.2.  Therefore,  from  simulations  we 

cannot  obtain  the  stabilization  of  the  variance  to  the  log- 
irradiance  variance  of  a  spherical  wave.  Nevertheless,  the  results  of  stochastic  simulations  are  in  complete 
agreement  with  the  results  of  the  weak  fluctuation  theory  and  the  ’’effective"  theory.  Numerical  simulations 

show  that  for  all  Q.  the  log-irradiance  variance  cr^^(O)  in  the  Andrews  model  exceeds  the  log-irradiance 

variance  a^/(0)  in  the  von  Karman  model. 

In  strong  turbulence  =  5*  10*^^  cm'^^^  the  behavior  of  the  log-irradiance  variance  changes 
dramatically  (Fig.  6.12).  In  the  simulations  the  range  of  variations  in  the  Rytov  variance  p^  coincide  with  that 
considered  in  Ref  46.  When  Q  <  Q/  the  variance  C7^^(0)  climbs  more  rapidly  than  tlie  Rytov  variance  Pq  and 
exceeds  Pq  in  the  range  0.3<0<0.7.  Then  the  growth  of  slows  down  and  lies  below  Pq  .  The  value  of 
O,  for  which  ^l^j=Pl  is  close  to  the  parameter  corresponding  to  the  effective  beam  waist: 

Recall  that  due  to  the  shift  of  the  effective  beam  waist.  Simulation  results  are  close  to  tlie  predictions 

of  the  "effective"  theory  for  However,  in  contrast  to  these  predictions,  we  did  not  observe  the 

minimum  in  the  vicinity  of  klf.  The  further  analysis  is  needed  to  understand  this  discrepancy. 

For  the  collimated  beam,  propagating  in  weak  turbulence,  the  dependencies  of  the  log-irradiance  and 

Rytov  variances  on  zlie  quite  closely  (Fig.  6.13).  In  strong  turbulence  {pl>\)  stabilizes  to  a  value  1.8. 

Before  the  saturation  at  the  diffractive  beam  edge  is  noticeably  larger  than  o-|^^^(0),  while  in  the  region 

of  stabilization  the  difference  is  not  so  pronounced. 


Fig. 6.11,  The  dependence  of  the  log-irradiance 

variance  on  zl ka]  at  the 

diffractive  edge  and  at  the  center  of  the 
focused  beam,  propagating  in  weak  turbulence 

Cj'=5-10-‘"  cm'''\ 
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Fig.6.12.  The  dependence  of  the  log-irradiance  variance  on  z!  ka]  at  the  diffractive  edge 
and  at  the  center  of  the  focused  beann,  propagating  in  strong  turbulence  =  5- 
Fig. 6.13.  The  dependence  of  the  log-irradiance  variance  on  Q.-  z!  kal  at  the  diffractive  edge 
and  at  the  center  of  the  collimated  beam,  propagating  In  strong  turbulence 


The  performed  simulations  allow  to  make  the  following  conclusions: 

1.  In  weak  turbulence  the  irradiance  fluctuations  are  close  to  statistically  homogeneous.  The  log- 

irradiance  variance  at  the  diffractive  beam  edge  exceeds  the  variance  at  the  beam  center  by 

10-20%. 

2.  In  the  vicinity  of  the  effective  beam  waist  the  irradiance  fluctuations  in  the  beam  section  become  statistically 

inhomogeneous.  In  weak  turbulence  the  log-irradiance  variance  at  the  diffractive  beam  edge  exceeds 

the  variance  cr^^(O)  at  tlie  beam  center  by  several  orders  of  magnitude.  In  strong  turbulence  the 

difference  between  a^/0)  and  is  not  large, 

3.  Beyond  the  beam  waist  the  irradiance  fluctuations  are  again  close  to  statistically  homogeneous. 

These  conclusions  to  some  extent  coincide  with  the  predictions  of  phase  approximation  in  Huygens- 
Kirchhoffs  method^  and  the  "effective"  theoiy^. 

For  the  collimated  beam  the  results  of  Monte  Carlo  simulations  are  in  complete  agreement  with 
theoretical  predictions. 
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Part  II 

Optical  Simulation  of  Laser  Beam  Propagation  through  the  turbulent 

Atmosphere 


1.  Introduction 


A  laser  beam  in  atmosphere  is  distorted  due  to  refractive  index 
inliomogenuties.  Statistics  of  phase  distortions  of  plane  or  spherical  wave 
propagating  through  turbulent  medium  is  well  investigated ^  In  practice  we 
often  have  to  deal  with  a  more  complicate  situation  -  imaging  of  an  extended 
object  illuminated  by  a  laser  beam  in  an  optical  ranging  system.  In  such  cases 
computer  or  optical  simulation  of  laser  beam  propagation  through  distorting 
medium  is  helpful^. 

Let  us  consider  a  laser  beam  that  is  back  scattered  so  that  angle  0 
(Fig.l)  is  small  enough.  The  total  phase  deviation  induced  by  atmosphere  may 
be  represented  as  a  sum  of  forward  and  back  path  terms,  so 
L(p  = 

<  A(p^  >=<  >  +  <  Aq>l  >  4-2  <  A(p^A<P2  >=2o^(\+K^{6))  (1) 

0^  =<  A{p}  >=<  A(pl  >, 
here  , 

K^($)  =<  A(p^A(p2  >  Icr 

We  can  see  from  (1)  that  for  small  values  of  the  angle  6,  when  phase 
distortions  along  the  beam  path  1  -  R  and  2  -  R  are  almost  the  same,  K^{0)^  i 
and  <  A(p'^  >  is  twice  greater  than  that  for  large  angles,  when  K^{0)  « 1.  This 
phenomenon  is  known  as  amplification  of  phase  fluctuations  in  the  back 
scattered  beam  and  is  one  of  so-called  double-pass  effects.  The  theory  of  these 

effects^,  well  developed  for  such  ideal  scattering  objects  as  a  point  reflector  or 
a  plane  mirror,  becomes  more  and  more  sophisticated  for  complex  objects  and 
needs  computer  simulations  or  optical  modelling. 

Some  related  phenomena  appear  in  imaging  of  an  extended  object^.  Let  us 
consider  two  rays  coming  from  different  points  of  an  object  to  a  single  point  of 
receiving  aperture:  I-M  and  2-M  in  Fig.2.  Turbulence  induced  phase 
difference  between  these  two  rays,  important  for  the  analysis  of  so-called 
anisoplanatic  effects,  is  A(p  -  Acp^  -  A(p2  and  its  mean  square  value  may  be 
written  as 

<  Acp^  >-<  A(p]  >  +  <  A(p\  >  -2  <  A(p^A(p2  >=  2o^(l  (2) 

Correlation  function  in  (1)  and  K^iO)  in  the  last  formula  are 

the  same  for  a  statistically  homogeneous  layer,  but  differ  if  turbulence 
parameters  depend  on  the  coordinate  z  along  the  optical  path.  These  simple 
e.xamples  show  close  ties  between  double  pass  and  anisoplanatic  phase  effects. 
In  general  cases,  such  as  imaging  of  an  extended  object  by  a  finite  aperture 
optical  system  or  scattering  by  a  non  point  object,  phase  effects  are  more 
complicated  and  may  be  described  by  means  of  correlation  or,  in  more  general 

case,  mutual  structure  function^  (see  Fig.3); 

=<  >  (3) 

Both  K^{0)  and  provided  they  exist,  may  be  expressed  through 

Now  one  can  conclude  that  for  simulation  of  double  pass  or 
anisoplanatic  effects  we  should  make  a  model,  that  enable  us  to  obtain  the 
mutual  structure  function  with  desired  accuracy.  Note  that  such  a 


Fig.1 
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model  does  not  implement  amplitude  effects,  but  only  phase  ones.  In  many  cases  this  approach  is  sufficient  for 
adaptive  optics  systems  design  and  analyses  of  their  efficiency.  In  other  case  a  more  complicated  model  with  a 
number  of  phase  screens  should  be  used. 

2.  Theoretical  background 


It  is  known^  that  in  a  homogeneous  turbulent  layer  with  Kolmogorov  power  spectrum  mean  square 
phase  difference  between  two  parallel  rays  is 

<A(p^  >=2.92CyLr^^^  W 

where  is  the  refractive  index  structure  constant,  k  -  the  wave  number,  L  -  the  path  length  and  r  -  the  rays 
separation  distance. 

If  the  rays  are  not  exactly  parallel,  and  the  angle  0  between  them  is  small  enough  one  can  apply 
expression  (4)  to  each  thin  layer  (see  Fig.3)  and  summarise  their  contributions  along  the  path: 

<<p^  >=  2.92/c'  J  Cl (z)lr{z)f^ dz.  (5) 


In  the  case  considered  above 


r(z)  =  + 


If  does  not  depend  on  z  the  integral  in  (5)  can  be  calculated  both  anal>lically  (by  a  rather 
sophisticated  expression)  or  by  a  computer.  For  non-Kolmogorov  spectra  the  analysis  is  more  complicated  and 
a  different  technique  is  used^^"^.  The  integral  representation  (5)  of  the  mutual  structure  function  D  can  be 
easily  applied  to  compare  it  with  different  models  of  turbulent  layer.For  computer  and  optical  simulation  of 

imaging  system  performance  a  number  of  phase  screen  models  has  been  used^'^^.  These  models  were 
designed  to  simulate  single  pass  optical  wave  propagation  and  their  possibility  for  double  pass  and 
anisoplanatic  effects  simulation  was  not  examined  yet.  That  can  be  done  calculating  the  errors  in 
representation  of  mutual  structure  function  by  a  chosen  model.  Anohter  way  is  to  construct  a  model  of 
turbulent  layer  in  such  a  way,  that  it  should  reproduce  properly  correlation  between  two  rays,  propagating  at  a 
small  enough  distance. 

—  Let  us  consider  a  model  consisting  of  N  phase  screens  placed 

^  ^  ^  between  planes  A  and  B  at  arbitrary  positions  (Fig.4).  Each  phase  screen  is 

1  ^  /  described  by  its  structure  function  Dj(r)  =  Ajr^'^^.  Suggesting  these  phase 

1  /  screens  are  statistically  independent,  and  ApA  one  can  write  for  mean 

I  /  square  phase  difference  between  beams  l-T  and  2-2': 

2.  <p^>='L«pI>  =  a-zH^'‘'  _  /6) 

”■1 . /  The  formula  (6)  is  a  discrete  approximation  for  integral  (5),  and  the  problem 

‘■‘I . T  now  is  to  choose  properly  the  number  of  the  screens  and  their  coordinates  to 

1  /  provide  sufficient  accuracy.  In  more  general  case  screen  weights  A^  may  be 

I  /  not  equal. 

-|  1  12  consider  some  physics  principles  for  optimisation  of  the  set 

B  - ^ -  0  Qf  phase  screens: 

Fig.4  1.  The  distortions  of  plane  wave  should  be  represented  exactly.  For  this 

purpose  we  can  use  only  one  phase  screen,  or  a  number  of  arbitrary  located 
screens.  The  sum  of  their  intensities  Aj  should  be  equal  to  the  total  coefficient  for  the  plane  wave: 

X;4.  =  2.92.C„VL  =  ^  (7) 

The  relation  (7)  does  not  define  different  values  Aj  and  coordinates  Zj,  but  it  should  be  valid  for  any  phase 
screens  representation. 

2.  The  model  should  represent  exactly  not  only  the  plane  wave  distortions,  but  also  distortions  in  spherical 
waves  emitted  both  in  planes  A  or  B.  As  the  configuration  is  symmetrical  we  may  try  to  design  a  system  of  two 
identical  phase  screens,  located  at  equal  distances  from  planes  A  and  B. 

Using  (7)  we  find  Ai=A2=Ao/2.  From  (6)  we  can  write 

^[[hef^^+[(L  -  h)ef^)  = 
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The  coefficient  3/8  in  the  right  part  of  the  last  equation  appears  from  the  well  known  result  for  phase 
fluctuations  in  a  spherical  wave^. 

This  gives  us  an  equation  +  (1  -  =  3/4,  a=h/L.  (8) 

The  numerical  solution  of  (8)  gives  a=0.21.  A  regular  technique  for  phase  screen  model  design  can  be 
based  on  the  Gaussian  quadrature  formula that  enables  us  to  calculate  best  coordinates  Zj  and  screen 
weights  for  a  given  number  of  screens.  For  two  screens  the  result  obtained  by  the  last  technique  does  not  differ 
significantly  from  that  obtained  above. 


3.  Two  phase  screen  method 


Let  us  consider  a  model  of  turbulent  layer  consisting  of  two  phase 
screens  placed  as  shown  in  fig. 5.  For  mutual  staicture  function  one  can 
write 

n)  =  ^  {[(1  - 

Applying  a=0.21  we  obtain  an  analytical  approximation  for  D^p  : 

0«,{r„rj)  =  A  j[o.79F,  +  0.21fjf '  +  [0.21?,  +  0.79?j]’^^}  (10) 

here  Aq  =  2,92  L,  (see  the  coefficient  before  r^^^  in  the  formula  (4) 

for  mean  square  phase  fluctuations  in  a  plane  wave). 

The  last  formula  is  useful  both  for  analytical  consideration  and 
computer  simulation  of  anisoplanatic  effect  in  statistically  homogeneous 

layer.  If  dependence  on  tlie  coordinate  z  along  the  path  is  significant. 


Fig.6 

Approximation  ererors  of  calculations:  A  -  rj  1^2 ,  B  - 


The  single  screen  model  with  double  beam  pass 
and  reflecting  random  phase  screen. 


the  coordinates  of  the  phase  screens  should  be  calculated 
using  Gauss  quadrature  technique. 

We  have  compared  the  analytical  approximation 
(10)  with  result  obtained  by  direct  computation  of  £)^(?j,^) 

via  integral  representation  (5).  The  difference  5d  between 
these  two  solutions  is  shown  in  Fig.6a  (for  q|  and  in 
Fig  .6b  (for  /jJU^).  We  can  see  that  in  the  worst  case 
approximation  error  does  not  exceed  8%.  In  other  cases  it  is 
much  less.  The  function  D^{T\,r^)  is  most  frequently  used 
in  integral  expressions,  so  the  accuracy  obtained  is 
sufficient  for  many  applications.  Two  screen  model  is  useful 
particularly  for  optical  modelling,  because  it  is  very  difficult 
to  construct  such  a  model  with  more  than  two  screens. 
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For  simulation  of  short  scale  distortions  a  single  screen  model  with  double  beam  pass  can  be  used.  In 
Fig.7  a  possible  optical  scheme  with  reflecting  random  phase  screen  is  shown. 

4.  Optical  simulation  technique  for  anisoplanatic  and  double  pass  effects 
investigation 

A  phase  screen  method  can  be  applied  to  design  an  optical  setup  for 
turbulence  induced  phase  effects  investigation.  To  simulate  "frozen"  turbulence 

computer  designed  phase  screens  may  be  used^^.  For  optical  “on-the-bench” 
turbulence  simulation  light  controlled  spatial  phase  modulators,  such  as  liquid 

crystal  light  valve  (LCLV),  are  of  great  interest^* These  devices  are  especially 
useful  to  produce  phase  screens  with  relatively  short  scale  small  amplitude  phase 
distortions.  As  we  can  see  from  previous  consideration  this  ability  is  of  great 
importance  in  anisoplanatic  and  double  pass  effects  studying.  As  for  large  scale 
distortions,  especially  low  order  Zernike  modes,  they  can  be  successfully 

simulated  by  computer  controlled  flexible  mirrors^’^. 

A  schematic  of  the  possible  optical  arrangements  is  shown  in  Fig.  8.  A 

main  difference  between  this  optical  scheme  and  that  proposed  before^  is  that  the 
light  beam  generated  by  laser  source  Ls  is  twice  reflected  by  LCLV  surface  in 
accordance  with  conclusions  of  the  theoretical  analysis  given  above.  Another 
Fig. 8  difference  is  that  to  produce  a  random  phase  modulation  a  laser  induced  speckle 

Schematic  of  the  possible  field  is  supposed  to  control  LCLV.  A  flexible  mirror  4  can  be  used,  if  needed,  to 
experimental  induce  low  order  phase  distortions.  After  reflection  from  this  mirror  the  laser 

arrangements  beam  reaches  the  input  plane  5  of  an  imaging  system,  designed  for  the  performed 

experiment.  An  image  of  a  test  object  2  or  intensity  distribution  in  output  laser 
beam  is  registered  by  a  CCD  camera  and  processed  by  a  computer.  The  key  element  of  our  experimental  setup 
is  a  reflective  LCLV  3,  that  transforms  the  intensity  of  the  speckle  field  into  random  phase  modulation  of  the 
redout  beam. 

The  general  configuration  of  the  device  is 
shown  in  Fig. 9.  The  ac  light  valve  consists  of  a 
number  of  thin  film  layers  sand^viched  between  two 
glass  substrates.  A  low  voltage  (5-20  Vj-n^s  )  audio 
frequency  power  supply  is  connected  to  the  two 
outer,  thin  film  indium-tin-oxide  (ITO)  transparent 
electrodes  and  thus  across  the  entire  thin  film 
sandwich.  The  50-micron  gallium  arsenide  (AsGa) 
semiconductor  is  used  as  photosensitive  layer.  The 
dielectric  mirror  and  cadmium  telluride  light 
blocking  layer  isolate  photoconductor  from  the 
readout  light  beam.  The  optical  isolation  between 

readout  light  and  photoconductor  exceeds  10^. 
This  is  one  of  the  major  design  features  of  the  ac 
light  valve,  as  it  enables  simultaneous  writing  and 
reading  of  the  device  without  regard  to  the  spectral 
composition  of  the  two  light  beams.  Furthermore, 
the  dielectric  mirror  prevents  the  flow  of  the  dc 
liquid  crystal  mirror  current  through  the  liquid  crystal  and  this  enhances 

+  CdTc  the  lifetime  of  the  device.  The  chemically  inert 

pjg  9  insulating  layers  Si02  bounds  the  liquid  crystal 

Liquid  crystal  spatial  phase  modulator  layer  at  both  interfaces. 

The  liquid  crystal  that  is  used  in  this 
device  is  a  biphenyl  nematic  material.  The  thickness  of  that  layer  is  10  microns.  The  optical  birefringence 

effect  of  nematic  liquid  crystal is  used  in  this  device.  The  liquid  crystal  molecules  at  the  electrodes  are 
aligned  with  their  long  axes  (director)  parallel  to  the  electrode  surfaces  (homogeneous  alignment).  In  addition, 
they  are  aligned  to  lie  parallel  to  each  other  along  the  preferred  direction  that  is  fabricated  into  the  device.  In 
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Phase  Delay  (rad) 


The  sensitivity  curve  of  the  LCLV 

order  to  obtain  pure  phase  modulation,  the  directions  of  liquid  crystal  alignment  are  the  same  at  the  both 

interfaces  (instead  twisted  nematic  used  in^^).  For  the  same  purposes  the  director  should  be  oriented  parallel 
with  the  readout  beam  polarization. 

Without  input  illumination  the  dynamical  resistance  of  the  photoconductor  layer  under  applied  ac 
voltage  is  much  higher  then  the  resistance  of  liquid  crystal.  Consequently,  applied  voltage  almost  completely 
falls  on  the  photoconductor.  Let  us  consider  what  happens  in  a  illuminated  element.  In  the  idealized  case,  the 
incident  photons  introduce  a  leakage  resistance  across  the  photoconductor  and  a  part  of  applied  voltage  acts  on 
the  liquid  ciy^stal  layer.  This  voltage  causes  the  S-type  deformation  in  the  liquid  crystal  layer,  and  as  a  result 
the  local  refraction  index  variation. 

The  response  time  of  the  device  is  20-25  ms  for  phase  modulation  equal  tc,  X=0.63|am.  The  spatial 
resolution  is  about  10  line  pairs/mm  at  the  50%  MTF.  The  sensitivity  curve  for  LCLV  is  presented  on  the 
Fig.  10.  We  can  see  that  in  the  interval  (0-27i)  phase  modulation  is  almost  linear  via  input  intensity. 

5.  Random  optical  field  generation. 


For  generation  of  time-independant  or 
dynamic  random  optical  field  speckle  technique 
may  be  used,  see  fig.  11.  Optical  scheme  includes 
lenses  LI,  L2  and  diaphragm  D  to  focus  laser  beam 
at  the  receiving  edge  of  a  fiber  bundle  FB.  A  mask 
M  can  be  placed  in  front  of  the  focusing  system  in 
order  to  control  the  light  intensity  distribution  at  the 
FB  edge.  The  receiving  edge  of  FB  is  placed  on 
adjustable  translating  mount,  that  allows  to  change 
the  size  of  the  light  spot  entering  the  FB  edge.  The 
second  edge  of  the  bundle  may  be  considered  as  a 
set  of  elementary  coherent  light  sources.  The  phases 
of  these  sources  are  randomly  distributed  in  (0,2;t) 
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Fig.12 

Optica!  speckles  patterns 

interval  and  their  amplitudes  corresponds  to  the  intensity  distribution  at  the  input  edge  of  FB.  The  speckle  field 
at  the  screen  S  is  registered  by  CCD  camera.  Patterns  of  random  light  intensity  distribution  for  a  number  of  FB 
input  edge  positions  are  shown  in  fig.  12. 

The  electronics  used  in  experiments  allows  us  to  obtain  digital  intensity  distribution  images  with 
512x512  resolution  and  256  gray  scale  coding.  The  total  field  of  view  on  pictures  shown  in  fig.  12  is  about  1.5 
cm,  so  we  can  easily  produce  speckle  fields  with  the  typical  avarage  spot  size  from  less  then  1  mm  to  about 
than  1  cm. 

First  step  to  calculate  cross-correlation  function  of  the  speckle  field  produced  is  smoothing  of  the  the 
intensity  array  ly.  In  order  to  surpress  high  spatial  frequency  noise  of  the  CCD  camera  and  relay  electronics  we 

smooth  picters  by  a  simple  mask: 

(l  2 

±2  4 

‘Hi  2 

To  get  tlie  estimate  of  the  speckle  field  correlation  function  we  select  a  number  of  equidistant  rows 
(usually  16)  and  calculate  the  “horizontal”  correlation  functions  in  accordance  to  the  formula: 

olz  A 

where  by  we  denote  the  avarage  intensity  of  the  i-th  row.  Almost  the  same  formula  is  used  to  calculate 
“vertical”  correlation  functions  Kyj  (avaraging  along  the  colomns  of  the  picture).  The  larger  the  separation  A 
the  fewer  number  of  overlapped  points,  hence  the  higher  the  variance  of  the  estimate.  This  means  that  we  can 
only  make  reasonable  estimate  of  the  correlation  functions  at  separation  out  to  about  1/4  of  the  image  size,  so 
in  our  calculation  we  use  [0,128]  separation  range  for  A.  As  the  last  step  the  average  correlation  functions 
Ky,  k  are  calculated: 

^=^(4+4) 

16  /=!  16  ^=1  2 

Before  calculating  ‘general*  function  K,  the  “vertical”  correlation  function  Kyj  is  rescaled  to  the 
same  scale  as  Correlation  functions  for  the  images,  represented  in  fig.  12,  are  shown  in  fig.  13.  We  can  see 
that  k^  and  ky  are  very  similar,  so  the  speckle  field  produced  is  spatially  homogeneous  and  the  function  K 
may  be  considered  as  a  general  characteristics  of  the  field. 

It  is  a  well  known  fact  that  correlation  function  of  the  speckle  field  intensity  may  be  considered  as 
Fourier  transform  of  intensity  distribution  in  the  emitting  plane^^*^^^  Provided  the  input  plane  of  FB  is  the 
focal  plane  of  tlie  lens  system  L1,L2,  tlie  field  distribution  in  that  plane  is  proportional  to  Fourier  transfonn  of 
the  mask  M  transmittance  profile^^.  When  the  mask  is  binary,  that  is  one  part  of  it  transmits  light  completely 
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and  the  other  part  doesn't,  we  can  speak  just  about  its  shape  and  size.  Making  use  of  Fourier  transform 
properties  one  can  obtaine  the  next  formula  for  output  correlation  function: 


K{p)  -  c|  M{r  -\-  —  mp)  M{r  -  —  mp)d^r 


(11) 


where  m=iyf  is  the  scale  factor,  p  is  the  vector  of  a  point  in  tlie  output  plane,  and  r  is  the  vector  of  a  point  in 
the  mask  plane  and  c  -  normalizing  constant. 

While  deriving  (11)  it  was  assumed,  that  the  illuminating  beam  is  uniform  and  tlie  input  edge  of  FB  is 
located  exactly  in  the  focal  plane  of  the  L1-L2  lens  system.  A  small  shift  of  the  input  edge  of  FB  complitely 
changes  the  correlation  function,  and  unfortunately  tlie  expression  for  K(p)  in  this  case  becomes  much  more 
complicated. 

In  practice  we  can  obtain  desired  shape  of  the  generated  speckle  field  correlation  function  by  changing 
the  mask  shape  and  adjusting  the  bundle  input  edge  position.  Applying  this  random  light  field  to  the 
photosensitive  layer  of  LCLV  we  can  produce  optically  controlled  phase  distortions. 

Dynamical  phase  distortions  can  be  simulated  by  moving  the  fiber  bundle.  Small  displacement  of 
emitting  end  provides  a  moving  speckle  pattern,  while  bending  its  middle  part  causes  a  change  of  random  field 
realisations,  somewhat  like  optical  turbulence. 


6.  Experimental  results 
on  random  phase  beam 
generation 

Optical  setup  made  in 
Adaptive  Optics  laboratory  to 
introduce  the  short  scale  phase 
distortions  in  laser  beam  is  shown 
in  fig.  14.  The  light  beam  of  the 
laser  Lsj^  is  expanded  by  the 

lenses  Li  X-2-  After  reflection  from 
the  LCLV  readout  surface  the 
phase  modulated  light  is  diverted 
by  the  beam  splitter  BSi  into  the 
neasuring  arm,  Sound  frequency 
generator  SFG  served  as  ac  power 
supply  for  LCLV. 

Control  speckle  field  is 
produced  by  laser  Ls2,  lens  system 


52 


Ls  and  fiber  bundl  FB.  Part  of  the  control  speckle  beam  is  diverted  by  the  beam  splitter  BS2  to  the  screen  S2 
and  video  camera  CCD2. 

To  make  the  phase  distortions  visible  we  used  a  Foukot  knife  technique^'^,  implemented  by  lenes  L3,L4 
and  knife  edge  K.  To  produce  the  image  of  the  LCLV  surface  on  the  screen  SI  the  distances  LCLV  -  L3,  L3  - 
K,  K  -  L4  and  L4  -  Si  equel  to  double  focul  length  of  lenses  L3,L4. 

Some  examples  of  control  speckle  field  and  produced  phase  maps,  visualized  by  Foukot  knife 


technique’  are  shown  in  Fig,  15,  corresponding  correlation  functions  for  this  patterns  are  shown  in  fig.  16. 


Control  speckle  field  (A,  B)  and  produced  phase  patterns  (C,  D) 


We  can  see,  that  correlation  functions  calculated  for  obtained  phase  patterns  reveals  good  coincidence 
with  that  of  the  control  fields.  Some  increase  of  the  short  scale  speckles  weight  in  the  phase  correlation 


functions  is  due  to  Foukot  knife  technique  influence. 


7.  Conclusions 


Possibility  of  double-pass  and  anisoplanatic  effects  simulation  by  few  phase  screen  model  is  discussed. 
A  two  phase  screen  model  of  a  turbulent  layer  is  proposed  and  its  accuracy  in  representation  of  mutual 
structure  function  is  examined.  An  experimental  set-up  for  experimental  investigation  of  anisoplanatic  effects 
witli  a  dynamic  phase  spatial  modulator  is  designed.  A  method  of  generation  of  random  optical  field  with  wide 
range  correlation  length  is  proposed  and  experimentally  tested. 
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Partin 

Photorefractive  Crystal  (PRC)  Based  Technology  for  Distortion 

Mitigation  in  Laser  Beams 


1.  Mitigation  of  phase  distortions  by  means  of  one-way  mitigation  system  with 
formation  of  pumping  wave  from  the  distorted  signal 

LL  Introduction 

During  some  last  years,  the  problems  of  formation  of  high-quality  laser  beams  and  compensation  or 
correction  of  distortions  resulting  from  propagation  of  such  beams  through  inhomogenious  medium  are  of  great 
interest.  These  problems  are  very  significant  in  such  important  applications  as  imaging  through  turbulent 
atmosphere,  laser  location,  transmission  of  laser  beam  through  optical  fibers,  etc.  The  main  goal  of  this  paper  is 
to  present  our  last  results  in  simulation  of  so-called  one-way  compensation  system  based  on  four-photon 
processing  of  distorted  signal  in  nonlinear  medium. 

The  first  systems,  used  to  correct  the  phase  distortions,  were  based  on  so-called  active  techniques.  In 
such  a  system,  wave  front  of  signal  is  changed  by  means  of  tlexible  adaptive  mirrors  ‘  or  holograms  before  its 
propagation  through  inhomogenious  medium.  A  control  of  the  signal  is  performed  in  such  a  manner  to  obtain 
required  wave  front  of  the  signal  after  its  way  through  the  medium.  Evidently,  for  optical  signals  with  rather 
complicated  wave  fronts  or  for  inhomogenities  that  change  in  time,  practical  realization  of  such  kind  a  system  is 
very  difficult  or  impossible. 

Another  way  to  solve  the  problem  lies  in  application  of  so-called  phase  conjugation  technique.  In  such  a 
system,  wave  front  of  signal  after  its  propagation  through  inhomogenities  is  directed  to  special  unit  -  phase- 
conjugating  mirror  (PCM)  -  which  performs  operation  of  phase  conjugation  of  wave  front  of  the  signal.  After 
that,  conjugated  signal  is  transmitted  through  the  same  path  in  the  reverse  direction  and  the  initial  wave  front  of 
the  signal  is  reconstructed.  In  fact,  very  many  nonlinear  processes  can  be  used  to  realize  PCMs.  However,  the 
most  promising  units  of  such  kind  are  usually  based  on  degenerate  three-  (2®  -  ca  =  ®)  or  four-  (®  +  ®  -  ®  =  ®) 
photon  processes  in  nonlinear  medium  ^  \  Here  ®  is  frequency  of  the  signal  and  some  additional  waves  -  pumping 
beams.  These  processes  are  tresholdless  and  enable  one  to  conjugate  ultrashort  laser  pulses  with  duration  up  to 
10  -10  s  .  Quantum  efficiency  of  such  PCMs  depends  on  value  of  nonlinear  susceptibility,  length  of 
nonlinear  interaction,  intensities  of  pumping  beams  and  increases  with  increase  of  these  parameters 

The  most  important  characteristic  of  any  mitigation  system  is  quality  of  distortion  compensation.  This 
characteristic  can  be  estimated,  for  example,  through  so-called  overlapping  integral 
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{p,z)dp 
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Here  Uq  and  f/om  are  the  amplitudes  of  the  input  (initial)  signal  and  the  output  (reconstructed)  field,  p  and  z 
define  position  of  each  point  in  polar  system  of  coordinates.  Many  investigations  of  H  dependence  on  various 
characteristics  of  input  signal,  inhomogenious  medium,  and  PCM  were  performed  It  was  shown  that  the  most 
promising  nonlinear  materials  for  the  practical  realization  of  PCM  are  photo-refractive  crystals  (PRC), 
Unfortunately,  all  schemes  of  such  kind  are  of  so-called  double-way  character.  That  means  that  phase  distortions 
resulting  from  the  first  propagation  of  signal  through  inhomogenities  can  be  compensated  only  after  one  more 
path  of  the  conjugated  signal  through  the  same  inhomogenities.  Limitations  resulting  from  this  specific  feature  of 
such  schemes  are  most  evident  in  the  case  with  rather  long  path  of  the  signal  through  inhomogeneities  that 
change  in  time.  In  this  case,  the  problem  of  mitigation  of  distortions  can  not  be  solved  by  such  a  manner.  That  is 
why  so-called  one-way  compensation  systems  were  proposed  Unfortunately,  up  to  now  many  problems  of 
their  practical  realization  and  operation  have  not  been  investigated  and  solved. 

One-way  compensation  system  can  be  realized  by  different  ways.  For  example,  one  can  pass  signal  beam 
Us  with  useful  information  and  pumping  beam  Upi  with  plane  wave  front  through  the  same  inhomogenities 
(fig. 1. 1).  In  case  when  these  inhomogenities  are  of  "pure"  phase  character,  one  can  subtract  them  from  wave 
front  of  output  wave  by  means  of  four-beam  interaction  in  PRC 
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Fig  1.1  One-wave  mitigation  system  with  pumping  wave  Up^  formed  by  the  same  inhomogenities  from  the 
special  wave  Up. 


Here  exp(i(j>)  is  the  transmission  function  of  inhomogenities.  It  is  easy  to  see  that  if  the  second  pumping  wave  Uf2 
is  plane,  four-photon  mixing  of  these  beams  in  nonlinear  medium  gives  rise  to  output  signal  without  any 
distortions. 

Unfortunately,  a  scheme  of  such  kind  always  needs  in  special  pumping  wave  Upi  transmitted  through  the 
same  path  as  the  signal.  Obviously,  it  is  not  so  convenient  in  many  applications.  However,  in  case  with  different 
spatial  scales  of  useful  signal  and  inhomogenities,  it  is  possible  to  use  another  way.  Processing  some  part  of 
distorted  signal  one  can  form  pumping  beam  with  needed  shape  of  wave  front.  To  do  this,  the  signal 
information  must  be  previously  removed  from  this  part  of  distorted  signal  by  special  spatial  filter  (fig.  1.2)  while 
information  about  all  inhomogenities  must  be  retain.  After  that,  the  distorted  signal,  the  pumping  wave  prepared 
by  this  way,  and  one  more  pumping  wave  with  plane  wave  front  must  be  mixed  in  PRC. 


Fig  1  2  One-wave  mitigation  system  with  formation  of  pumping  Upi  by  spatial  filtration  of  distorted  signal 
Us. 

Two  main  factors  can  influence  on  quality  of  compensation  of  phase  distortions.  The  first  factor  is 
connected  with  accuracy  of  phase  conjugation  and  is  inherent  to  all  versions  of  such  a  system.  This  procedure  is 
not  ideal  and  its  accuracy  is  determined  by  PCM  transfer  function.  Efficiency  of  this  nonlinear  process  depends  on 
pumping  intensities,  the  nonlinearity,  the  length  of  nonlinear  interaction,  etc.  At  the  same  time,  self-action  and 
stimulated  scattering  as  well  as  limited  angular  aperture  of  PCM  will  decrease  quality  of  conjugation  of  input 
signal.  The  second  factor  is  inherent  to  the  version  of  one-way  mitigation  system  with  spatial  filtration  of  the 
distorted  signal.  Removing  of  a  part  of  signal  information  and  retaining  of  a  part  of  information  about 
inhomogenities  result  in  decrease  of  quality  of  distortion  compensation.  The  best  quality  must  correspond  to  the 
case  when  shape  of  wave  front  of  formed  pumping  beam  coincides  with  shape  of  wave  front  of  special  plane 
pumping  beam  distorted  by  inhomogenities.  However,  because  spatial  spectra  of  transmitted  information  and 
inhomogenities  overlap,  it  is  not  possible  to  form  such  ideal  pumping. 

A  possibility  of  distortion  mitigation  by  such  a  scheme  can  be  illustrated  by  a  simple  example  relating  to 
the  case  with  discrete  spectra  of  signal  and  inhomogenities.  Let  us  suppose  that  transmitted  informative  signal 
is  given  in  the  form  of  one-dimensional  harmonic  phase  grating 


U,  =  Uq  exp{/acos(ax:))  =  UoZ/'”y„(a)exp(//?7av).  ^13^ 

Here  Uq  is  the  amplitude  of  the  signal,  a  and  x  represent  the  depth  and  spatial  frequency  of  the  useful  phase 
modulation.  Let  us  suppose  that  phase  distortions  can  be  described  as  similar  phase  grating  with  the  amplitude  b 
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and  spatial  frequency  p.  That  means  that  after  propagation  through  the  inhomogenious  medium,  the  useful  signal 
multiply  by  factor  exp(ibcos(Px))  and  can  be  represented  as 


U[  -  S/""'”y^(a)y„(Z?)exp(/(wa  +  /?>^)x). 

Therefore,  its  spatial  spectrum  which  is  shown  in  fig.  1.3  can  be  written  in  the  form 

J ^{a)j ^{b)5[ka  +  /?/?+/<:). 

k,p=-^ 


Fig.1 .3  Spatial  spectrum  of  the  signal  Us  distorted  by  harmonic  phase  grating. 

Ideal  pumping  wave  front  (see  above)  must  include  all  spectral  components  with  k^k^.  Here  Jico(o)  =  max  Ji^  (a), 
k  =  0,±],  ±2,,..  .  Corresponding  spectrum  of  pumping  is  shown  in  fig.  1.4.  Quality  of  mitigation  will  decrease  by 
two  different  reasons.  The  first  reason  results  from  presence  of  some  pumping  spectral  components  with  k  ^  ko 
and  the  second  one  follows  from  absence  of  some  components  with  k=kQ. 


Fig.1 .4  Spatial  spectrum  of  the  pumping  wave  Upi  formed  by  spatial  filtration  of  the  distorted  signal  UV 

This  situation  is  typical  for  any  real  conditions  and  must  be  simulated  in  calculations.  Another  problem  results 
from  random  behavior  of  real  inhomogenities  which  are  of  interest.  That  means  that  stochastic  description  must 
be  used  to  optimize  one-way  mitigation  system. 

1,2.  Description  of  the  model 

Our  theoretical  model  of  proposed  one-way  mitigation  system  is  based  on  a  step-by-step  analytical 
description  of  the  input  image  transformation.  That  means  that  in  our  model,  we  consequently  consider  the 
transformation  of  the  input  image  after  all  parts  of  its  total  way  through  the  homogeneous  region  of  atmosphere, 
the  turbulent  region,  the  filter  of  spatial  frequencies  for  the  channel  of  the  pumping  wave  formation,  the  PCM, 
etc.  (fig.  1.5).  On  each  such  a  step  we  use  the  paraxial  approximation. 

For  example,  after  the  way  through  the  homogeneous  region  of  atmosphere  with  the  length  Az,  spatial 
spectrum  of  the  signal  U  (k)  transforms  as 


f/(/c,z  + Az)  =  U[k,z) 


1 

IK^ 

~2k 


\ 

Z\z 

/ 


(1.6) 
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Here  k  and  /care  the  magnitude  and  the  transverse  projection  of  the  wave  vector  of  the  signal.  We  consider  only 
two-dimensional  case,  but  the  model  can  be  easily  generalized  to  three-dimensional  case  ifx  and>^  coordinates  are 
independent. 
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Fig. 1.5.  Block-diagram  of  one-way  mitigation  system:  Lf  -  lenses  of  spatial  filter;  Li,2  -  lenses  focusing  the 
signal  and  pumping  into  PCM;  L  -  lens  focusing  the  restored  signal  on  the  plane  ze. 

We  simulate  turbulent  region  of  atmosphere  by  a  set  of  u  thin  random  phase  screens  with  transmission 
functions  exp[/V//')].  Here  <l>^  is  the  random  phase  shift  of  the  signal  after  its  path  through  the  /7-th  screen. 
Therefore,  we  describe  the  input  signal  transformation  by  each  the  screen  as 

U'{r,z„)  = 

Here  is  the  longitudinal  coordinate  of  the  /7-th  screen.  Expression  (1.7)  is  equivalent  to  transformation  of 
spatial  spectrum  of  the  input  signal  which  is  given  by 


^7'(/c,zJ=  \  u[k\z yi\{K  -  K’)dK 


(1.8) 


Here  TJ^k-k^  is  the  Fourier  image  of  the  transmission  function  exp[i(l>n(r)].  We  consider  the  way  of  the  signal 
between  each  pair  of  phase  screens  as  homogenious  region  of  atmosphere  and  described  the  input  signal 
transformation  on  these  parts  of  the  total  way  by  expression  which  is  similar  to  (1.6).  Such  the  model  is 
reasonable  if  distance  between  the  screens  {A)  are  more  than  the  correlation  length  in  z  direction.  This 
assumption  enables  us  to  consider  the  screens  as  independent. 

We  consider  the  filter  of  spatial  frequencies  to  be  composed  of  two  lenses  with  the  same  focal  length  ff 
and  a  pin-hole  placed  in  the  middle  of  the  distance  2ff  between  these  lenses.  The  filter  transforms  spatial 
spectrum  of  the  signal  from  the  input  focal  plane  z.f  of  the  first  lens  to  the  output  focal  plane  z+f  of  the  second  lens 
in  accordance  with  expression 


/  \  h  !  \  (  J'fA 


(1.9) 


Here  /f  is  transmission  function  of  the  pin-hole.  Further,  the  planes  z.„  and  z=n  are  considered  as  a  plane  Z2 
(fig.  1 .5  ).  It  means  optical  ways  of  the  signal  and  pumping  are  the  same. 

The  PCM  is  assumed  to  be  an  ideal  four-photon  converter  of  the  input  fields  in  relation  to  its  spatial 
resolution  and  the  speed  of  response.  The  second  pumping  wave  is  assumed  to  be  plane.  In  this  case,  the  PCM 
performs  the  convolution  of  spatial  spectra  of  the  input  signal  and  formed  pumping  Up 


U'{k,z)  =  k]u^{k\z)uI{k’- K,z)dK’ 


(1.10) 


Here  K  represents  energy  efficiency  of  the  PCM. 


59 


An  additional  optical  system  combined  of  the  lenses  Lj  and  L2  with  the  same  focal  length/  and  the  lens 
L  with  the  focal  length  fi  is  used  to  focus  the  pumping  wave  and  the  distorted  image  into  PCM  and  to  form  the 
output  signal.  Expression  for  spatial  spectrum  of  the  output  (restored)  signal  can  be  obtained  in  the  form 

J  =  Cin. \.JdK^^dK^,,dK^,dKjK^,_ciKjK^clKdKU^  (  )y  ‘  ( ) 

h  (- ^3 )  -  Si  )  '^2  ( '^,-3  -  S2 )  'I'x  ( S2  -  Si )  '^’2”  ( 'f.,3  -  '^.2 )  exp(/«-) 


J  2  2  .  2 

pSiS  -S2^-S3 


A-S^(s-^i  -s)  + 


+U'c  -  S3  /  +  S3-,  '  4  A  '-s' 


3 Is  --1  -  A)  '-is  -S  Is  --3J 


J+(/C-/f)V2]  I 


Here  Uq  represents  spatial  spectrum  of  the  transmitted  image,  C  is  the  constant  which  depends  on  K,  k,  /  fi. 
Expression  (1.11)  relates  to  the  model  which  includes  only  two  random  phase  screens  but  any  number  of  such 
phase  screens  can  be  included  in  (1.11)  with  easy.  Furthermore,  we  suppose  that  (pi Jr)  corresponds  to  the 
isotropic  random  process  with  Gaussian  statistics.  We  suppose  the  distance  between  the  screens  are  more  than 
the  correlation  length  along  z  direction  (see  above).  This  assumption  enables  us  to  simplify  expression  (1.11)  by 
means  of  its  averaging  and  usage  of  relation 

(7;(  /c,)  /;( x-J  7;’(  x-,)  7;(  a'J)  =  5,(  /c,).s;(  x-,  -  k^^-k,). 


Finally,  we  find  expression  for  the  mean  value  of  the  output  signal  in  the  form 


(f/)= 


C\  .  .  J  (JKCh]Jl]jK^jK 

^''l(^l)‘^'2('72)v(-('fo  +  Vx  +  ^2)) 
-  A^)exp(7/cr) 

i\ A 

7  “^7: - 


k\  2 


exp]  -Jk  -  kY  -  k'‘{z,  -z  J  -  Af"(z5  -  zj] 
z/c 


Here  p  is  the  scale  constant  of  the  lens  system  Li  2, 


A^-z,-z,+[]  ^  JZ3-Z,  Z3  +  A 

=^4-23  +  A(23-^i). 
A„2=^4-^}+  fi[h-^X-^)- 


To  analyze  (1.14),  we  express  spectral  density  S(k)  in  the  terms  of  the  correlation  function  Bn  of  the  phase  (j) 

5'(Ar)  =  1  exp(-jKr-a‘ +.6i2(r))^/r,  I 

a' =77,2(0).  (1.1 
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To  find  variance  of  the  output  signal  we  have  to  average  expression 


exp 


1k5[k^  -  +  K^\ , .  cxp[iK^R^  +^^3^  “^^4^4) 

isjj!,) + -  A=(«.)  *A=(K, -i>.hBAK  -  k,)-bAr,  -  «,)! 


Expressions  (1 . 1 1)  and  (1 .20)  enables  us  to  find  variance  of  the  output  signal 

D=<U  U  *>-<U^,><U*>.  (1-21) 

out  out  out  out 

To  calculate  the  output  field  characteristics  in  analytical  way,  we  suppose  the  function  B,2  (r)  to  be 
Gaussian  and  expand  it  into  a  series 


-'•2)  = 


1- 


(1,22) 


Here  R  is  the  correlation  radius.  In  this  approximation,  the  spectral  density  .S'  is  a  Gaussian,  too, 

f 

S{T])  =  ^f^R^  exp 


(1.23) 


According  to  (1.22)  expression  (1.20)  is  equal  to 

{  ^1  (  )  ^1  (  ^^3  )  ^1  (  ^4  ))  “ 

=  6\(/cJb(/C2  -  fc,)B(/Cj  -  /cJbI/c,  -  /cj. 

Transmission  function  of  the  pin-hole  is  supposed  to  be  a  Gaussian,  too, 
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=  exp 
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(1.25) 

(1.26) 


where  c  is  the  pin-hole  radius. 

We  describe  the  input  signal  in  the  form  of  two  shifted  Gaussian  spots  and  suppose  the  first  spot  is 
located  in  the  middle  of  the  field  of  view  and  the  second  spot  is  shifted  in  relation  to  the  first  one  in  distance  ^2 


=  Ay  exp 

^  2  ^.2  ^ 
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(1.27) 


Here  A 12  and  Ro],o2  are  their  amplitudes  and  radii.  Such  shape  of  a  signal  enable  us  to  investigate  possibility  of 
image  reconstruction  over  the  field  of  view. 

Results  (1.14),  (1.21),  (1.22),  (1.25)  and  (1.27)  enable  us  to  obtain  mean  value  and  variance  of  the 
output  signal  in  analytical  form.  Because  final  expressions  are  rather  unwieldy  to  be  written  here,  they  will  be 
represented  here  only  in  a  graphical  form  (fig.l.6-fig.l.l5). 

L3.  Results  of  calculations 

At  the  first  glance,  it  seems  there  are  too  many  free  parameters  to  optimize  our  system.  However,  it  is 
not  right  and  we  can  choose  the  values  of  some  part  of  them  using  a  very  simple  qualitative  reasoning.  Let  us 
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consider  the  mitigation  system  with  only  one  phase  screen  and  plane  pumping  wave  passing  through  the  same 
inhomogenieties  as  the  informative  signal.  In  this  case,  spatial  spectrum  of  “ideal”  output  signal  will  be  defined  by 

K,  ( ^4 ))  =  K  {-pk)B{-AKlk)  exp(/4  pK  ^  Ilk).  (1-28) 

Here  B  is  the  Fourier  transform  of.V  (1.18),  that  is  a  correlation  function  of  exp(i(l>n).  It  follows  from  expression 

(1.28)  the  case  of  ideal  reconstruction  can  be  realized  only  if  both  lenses  L12  are  focused  onto  the  phase  screen 
plane.  That  means  the  distortions  in  the  PCM  plane  are  of  “pure”  phase  character.  In  this  case,  .^(0)=^!  and  the 

factor  Qxpi^lA^^K^  can  be  corrected  by  symmetrical  adjustment  of  the  system 


(1.29) 


It  is  easy  to  see  that  there  is  an  similar  requirement  for  the  system  with  two  phase  screens  and  the  plane  of 
optimal  focusing  lies  between  the  screens.  Its  position  is  determined  by  relationship  between  correlation  radii  of 
phase  distortions  In  the  case  when  the  radii  are  the  same,  the  plane  of  optimal  focusing  is  placed  in  the  middle 
of  interval  between  the  screens.  In  fact,  our  model  parameters,  which  should  be  optimized,  are  only 
characteristics  of  informative  signal,  inhomogeneities,  spatial  filter,  and  position  of  focal  planes  of  lenses  L. 
Therefore,  the  system  optimization  means  determination  of  optimal  relationships  between  them. 

To  check  up  the  ability  of  our  system  to  correct  phase  distortions,  we  performed  numerical  calculations 
of  analytical  expressions  for  mean  value  and  variance  of  output  signal.  Firstly,  it  was  necessary  to  find  out 
allowable  relationships  between  spatial  scales  of  informative  signal,  inhomogenieties,  and  the  width  of  spatial 
filter  that  enables  us  to  reconstruct  the  image.  For  these  relationships,  spatial  distribution  of  mean  value  of  the 
output  signal  must  coincide  with  initial  distribution  of  the  informative  signal.  Solving  this  problem,  we  considered 
some  different  combinations  with  one  (Ai  or  A2  =  0)  or  two  input  beams  and  one  or  two  phase  screens.  The 
lenses  L  were  focused  onto  the  plane  of  the  phase  screen  or  onto  the  middle  plane  between  two  identical  phase 
screens.  The  case  of  one  phase  screen  was  simulated  by  the  case  of  two  phase  screens  by  putting  the  distance 
between  them  equal  zero.  So,  in  the  case  of  one  phase  screen,  effective  correlation  radius  was  equal  to  R^j 4^ . 

To  filter  the  signal  spectral  components  and  to  retain  the  spectral  components  of  the  inhomogenities, 
one  needs  to  satisfy  some  requirements  to  their  spatial  characteristics.  Therefore,  one  can  expect  the  mitigation  of 
distortions  in  the  range 


(1.30) 


The  left  part  of  inequality  (1.30)  corresponds  to  passing  of  components  of  spatial  spectrum  of  the  distortions 
through  the  spatial  filter.  In  the  case  of  its  violation,  distortions  of  the  output  field  will  be  exactly  of  phase 
character  because  the  formed  pumping  wave  is  almost  plane.  At  the  same  time,  the  random  phase  distortions  do 
not  almost  subtracted  and  the  variance  of  the  output  field  must  be  very  large.  That  means  that  normalized 
variance  of  the  output  field  is  a  good  criterion  of  efficiency  of  this  subtraction.  The  right  part  of  inequality  (1.30) 
corresponds  to  filtering  of  some  components  of  spatial  spectrum  of  the  informative  signal  by  the  spatial  filter.  At 
the  same  time,  random  phase  distortions  are  subtracted  and  the  variance  of  the  output  signal  must  be  small. 
In  the  case  of  its  violation,  mean  value  of  the  output  signal  will  be  distorted  because  amplitude  profile  of  the 
formed  pumping  wave  will  repeat  the  amplitude  profile  of  the  informative  signal  and  the  last  will  be  squared.  It  is 
very  important  that  these  nonlinear  distortions  can  not  be  corrected  by  any  routine  optical  system.  For  the  case 
with  Gaussian  beam,  such  distortions  result  to  change  of  its  radius.  Therefore,  full  width  of  half  maximum 
(FWHM)  of  mean  value  of  the  output  signal  is  a  good  criterion  to  estimate  nonlinear  distortions.  Our  calculation 
enables  us  to  determine  an  effective  operating  range  of  one-way  mitigation  system  that  provides  the  output  signal 
of  both  minimal  nonlinear  distortions  and  minimal  variance. 

Results  of  our  numerical  calculations  of  the  output  beam  FWHM  and  its  variance  D  versus  and  the 
width  of  the  spatial  filter  /q  (1.26),  are  shown  in  fig.  1.6  and  fig.  1.7.  The  variance  is  normalized  to  the  unit 
magnitude.  The  calculation  results  are  in  full  accordance  with  the  above  discussion  of  inequality  (1.30).  It  is  easy 
to  see  that  in  the  case  when  the  informative  components  of  the  input  signal  spectrum  are  not  filtered  (/^n^o^^I)> 
the  output  signal  is  really  proportional  to  the  squared  informative  signal  and  the  limit  value  of  its  FWHM  is  equal 
to  V2/?q.  The  same  value  can  be  obtained  in  analytical  way  from  (1.14)  for  another  case  when  Rn«R{).  However, 
in  the  case  when  only  the  right  part  of  (1.30)  is  violated  i.e.  /^  <  1  /  R^  and  <\l  R^,  the  variance  D  is  small 
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(see  fig.  1.7),  too,  due  to  rather  large  spatial  scale  of  inhomogeneties.  To  determine  the  effective  operating  range 
of  the  mitigation  system,  v/e  calculated  the  isolines  (fig.  1.8)  of  F  =  [(FWHM-V2/^J  /  (IRq-^JiRq)]  (solid  lines) 
and  normalized  D  (dashed  lines)  which  are  plotted  on  the  \g{RJRo)  -  IgihRo)  plane.  The  step  between 
consecutive  two  isolines  is  equal  to  0,16.  One  can  see  that  in  the  case  of  F  >  0.32,  the  variance  D  <  0.16. 
Therefore,  the  major  limitation  on  effective  operating  range  of  our  system  is  determined  by  distortions  of  the 
mean  value  of  the  output  signal. 

Fig.  1.6-1. 8  refer  to  the  case  of  distortions  of  pure  phase  character,  that  is,  to  the  case  of  precise 
focusing  of  the  lens  system  onto  the  plane  of  two  coincident  phase  screens.  For  extended  atmospheric  path,  such 
a  case  can  not  be  realized.  That  is  why  we  calculated  FWHM  and  D  for  the  case  with  two  identical  phase  screens 
separated  by  distance  A  simulating  the  length  of  atmospheric  path.  Under  conditions  of  optimal  focusing  of  the 
lens  system  onto  the  middle  plane  (see  above)  and  transmitted  beam  radius  1  cm,  results  of  our  calculations  are 
shown  in  fig.  1.9-1.12.  Allowable  value  of  A  has  been  found  to  be  about  10^- 10^' cm  that  is  long  enough  for 
the  most  part  of  practical  important  applications.  The  major  limitation  on  effective  operating  range  of  mitigation 
system  is  determined  as  before  by  distortions  of  the  mean  value  of  the  output  signal. 

Formation  of  spatial  spectrum  of  the  pumping  wave  is  necessary  attended  with  its  angular  divergence. 
This  specific  feature  of  the  mitigation  system  can  result  in  dependence  of  the  PCM  characteristics  on  transverse 
coordinates  That  is  why  we  calculated  FWHM  and  D  of  the  output  field  in  the  case  of  the  one  input  Gaussian 
beam  (Ai  =  0)  shifted  by  R2  in  transverse  direction.  Calculation  results  are  shown  in  fig.  1.13  and  fig.  1.14.  It  is 
easy  to  see  that  quality  of  mitigation  does  not  practically  depend  on  initial  position  R2  of  the  input  beam. 

Calculating  the  case  of  two  input  Gaussian  beams  (Aj  =  A2),  we  checked  a  possibility  to  mitigate 
distortions  of  more  complicated  transmitted  image.  The  input  beams  were  supposed  to  be  originally  resolved  by 
the  Rayleigh  criterion.  Averaged  distribution  of  the  output  field  is  shown  in  fig.  1.15.  In  this  case,  criterion  of 
small  D  is  good  to  the  same  region  and  the  output  beams  can  be  resolved  even  better  than  the  input  ones  due  to 
some  nonlinear  distortions  (see  above).  The  Rayleigh  criterion  can  be  violated  only  for  very  large  A  ~  10^-10^  cm, 
that  is,  for  very  long  atmospheric  way. 

In  a  real  PCM  with  the  PRC  length  1,  nonlinear  interaction  should  be  considered  taking  into  account 
phase  mismatch  S  of  wave-vectors  of  interaction  waves  along  z-axis.  The  spatial  spectra  of  the  output  field  is 
determined  by  the  expression 


« 

U'{/c,z)  =  K\  \  -  fc,z')cxp{fSz')dz'djc\  (131) 

-1/2^  ^  ^ 

where  in  paraxial  approximation  S  =  ~  One  should  integrate  expression  for  the  output  field  Uou, 

obtained  for  the  ideal  case  (1  =  0)  over  the  PRC  length  because  parameters  and  p  should  be  replaced  by  the 
following  way 

A +2h+*/>^)' 

This  is  an  exact  way  to  take  into  account  the  real  PCM  length.  However,  much  more  simple  way  of  estimation  of 
such  kind  lies  in  calculation  of  an  equivalent  problem.  One  can  calculate  FWHM  and  D  for  the  case  when  the 
focusing  lens  system  is  mismatched  for  a  distance  about  1.  A  typical  value  of  1  is  about  only  1  cm.  It  is  much  less 
than  diffraction  length  of  the  input  signal  beam  with  Rq^  \  cm  and  can  not  considerably  influence  on  the  obtained 
results. 


1.4.  Conclusions 

1.  The  model  of  one-way  mitigation  system,  based  on  nonlinear  interaction  of  distorted  informative  optical 
signal  and  pumping  wave  formed  from  a  part  of  the  distorted  signal  has  been  developed. 

2.  A  criterion  of  successful  mitigation  of  phase  distortions  have  been  formulated. 

3.  Full  width  of  half  maximum  of  mean  value  of  the  output  field  and  its  variance  have  been  calculated. 

4.  An  effective  operating  range  of  the  mitigation  system  has  been  determined  and  optimized. 

5.  It  has  been  shown  the  optimized  mitigation  system  is  good  for  many  practical  important  cases. 
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2.  Dynamic  distortion  suppression  by  phase-conjugation  systems  based  on 
inertial  photorefractive  nonlinearity 

2. 1.  Introduction 

During  some  last  years,  the  problem  of  phase  distortion  mitigation  by  means  of  systems  based  on 
photorefractive  crystals  (PRC)  is  of  great  interest  \  The  main  reason  of  such  interest  is  a  possibility  to  realize 
high*efficient  nonlinear  interaction  using  very  low  intensity  of  pumping  beams  I  On  the  other  hand,  the  PRC 
nonlinear  response  is  inertial  and  nonlocal  ^  These  specific  features  enable  one  to  develop  some  new  and  very 
promising  techniques  to  suppress  the  dynamic  distortions.  In  conventional  version  of  mitigation  system  , 
nonlinear  medium  is  used  only  as  a  phase-conjugating  (PC)  mirror.  Therefore,  to  compensate  the  phase 
distortions,  one  needs  to  use  the  second  pass  of  conjugated  signal  through  the  same  inhomogenities.  This  pass 
greatly  narrows  down  the  field  of  possible  applications  of  such  a  system.  So  called  “one-way”  mitigation  systems 
allow  one  to  except  the  second  pass  of  restored  beam  through  the  inhomogenities.  In  systems  of  such  kind,  wave 
front  of  the  first  beam  -  the  signal  -  is  modulated  by  both  the  useful  information  and  phase  distortions.  The  second 
one  -  the  pumping  -  is  phase  modulated  too,  but  the  shape  of  its  wave-front  results  only  from  the  phase 
distortions.  To  restore  the  initial  information,  one  must  simply  subtract  the  wave  front  of  pumping  from  the  signal 
one  by  means  of  nonlinear  interaction  in  PRC.  To  form  the  needed  shape  of  wave  front  of  pumping,  one  can  use 
an  additional  plane  wave  transmitted  through  the  same  path^.  Another  possible  way  to  solve  this  problem  lies  in 
spatial  filtration  of  the  distorted  signal  Implementation  of  inertial  PRC  gives  us  the  new  possibility  to  realize  a 
crucially  new  type  of  one-way  system.  Such  a  system  is  based  on  averaging  of  dynamic  distortions,  that  is,  the 
nonlinear  medium  is  used  as  a  time  averaging  filter.  Authors  of  reference  ^  were  the  first  who  suggested  and 
realized  such  technique  and  published  some  preliminary  experimental  results.  The  most  promising  field  of 
applications  of  such  system  relates  to  problems  connected  with  transmission  of  optical  information  through 
turbulent  atmosphere. 

The  main  goal  of  our  paper  is  to  present  our  last  results  in  computer  simulation  of  dynamic  distortion 
suppression  by  the  one-way  system  based  on  nonlinear  interaction  in  photorefractive  crystals. 

2.2.  The  model 

We  took  the  conventional  set  of  material  equations^  as  a  basis  of  our  model.  Using  this  set  of  equations, 
we  obtained  all  analytical  steady-state  solutions.  To  compute  dynamics  of  all  the  processes,  we  used  the 
simplified  ^  version  of  the  set 
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(2.1) 


dic  dl  ’ 

j  =  +//,)(/%  +  £„)  +  /ukj-^, 


r  1  ^ 

-— =  al~a  n  + — 


Here  g  anci  T  are  the  static  dielectric  constant  and  the  temperature  of  PRC,  e  and  j-i  are  the  charge  and  the 
mobility  of  free  carriers,  is  the  Boltsman  constant.  Solution  of  (2.1)  remains  all  specific  features  of  dynamics 
of  spatial  distributions  of  the  electric  field  and  of  the  charge  p  and  the  current  j  densities.  We  supposed  that 
the  parameters  a  and  representing  photogeneration  and  recombination  of  free  carriers  are  the  constants,  that 
is,  they  do  not  depend  on  the  dark  and  photoinduced  densities  of  free  carriers.  We  described  the  external 
static  electric  field  applied  to  PRC  by  Eq.  So,  we  took  into  account  both  drift  and  diffusion  components  ofy. 

The  set  (2.1)  and  the  conventional  wave  equation  for  the  amplitude^  of  light  field  with  wave  vector  k 


31  , 

2ik  —  =A,A  +  2k^  —  A 
3  '  n 


(2.2) 


gave  us  the  complete  self-consistent  problem.  Here,  n  and  3i  are  the  linear  and  nonlinear  parts  of  the  refractive 
index.  We  calculated  the  last  term  by  means  of  well-known  electrooptic  equations^  and  solution  of  (2.1).  To 
simulate  the  problem  numerically,  we  represented  PRC  as  a  set  of  thin  nonlinear  layers  placed  consequently  along 
z-axis.  We  supposed  that  inside  each  layer,  all  functions,  including  the  beam  intensity  /,  depend  only  on  the 
transverse  coordinate  jc.  We  took  into  account  all  changes  along  z-axis  by  means  of  computing  (2.2).  We  assumed 
that  both  boundary  and  initial  conditions  correspond  to  the  system  '’switching  on"  at  the  moment  t=0  and  that 
both  the  interacting  beams  are  symmetrically  directed  onto  the  first  layer.  For  each  moment  /,  we  computed  (2.2) 
using  the  fast  Fourier  transformation.  Then,  we  made  the  step  in  time  and  calculated  Sfi  as  a  result  of  each  layer 
evolution  to  the  new  steady-state  distribution  of  all  the  parameters  listed  above.  In  this  procedure,  we  used  all 
initial  conditions  from  the  previous  time  moment.  So,  we  solved  this  self-consistent  problem  by  step-by-step 
technique. 

To  simulate  four-beam  interaction,  we  modified  our  numerical  scheme.  We  took  into  account  only 
transmitting  gratings,  that  is,  we  simulated  the  case  with  non-coherent  or  orthogonal  polarized  pumping  beams. 
This  assumption  enabled  us  to  speed  up  the  calculation  thanks  to  decrease  of  the  number  of  required  nonlinear 
layers.  The  main  specific  feature  of  our  numerical  scheme  for  this  geometry  was  the  following.  For  the  stage  of 
propagation,  using  inertia  of  PRC,  we  consequently  simulated  the  pass  of  all  the  beams  along  z-axis  in  direct  and 
reverse  directions.  Then  we  computed  the  nonlinear  part  <^7.  We  thoroughly  checked  accuracy  of  our  numerical 
scheme  by  changing  the  value  of  steps  of  the  grids  in  time  and  space. 

We  described  the  phase  distortions  arising  from  the  atmospheric  turbulence  as  some  additional  phase 
modulation  of  wave  front  of  the  signal  on  PRC’s  input  plane  (subscript  “in”)  according  to 

S<p‘''  =  /l,sin(/:^.)r)+ 4.sin(^  X -Q/).  (2.3) 

Here,  the  first  item  with  the  amplitude  A^  at  the  spatial  frequency  k^  represents  the  useful  information, 
and  the  second  one  with  A„  and  k^  describes  the  dynamic  phase  distortions  at  the  atmospheric  turbulence 
frequency  f2.  We  estimated  the  mitigation  quality  as  the  root-mean-square  deviation  of  the  restored  beam's  phase 
wave  front  (subscript  "out")  fro mThe  initial  information 

/?“'  =  -  A^s\n{kxfy. 

The  amplitude  distortions  were  simulated  as  an  additional  amplitude  noise,  that  is,  for  the  optical  field  in 
PRC’s  input  plane,  we  used  the  expression 

E'"  =  E"*™'+E„sin(k„x).  (2.5) 
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Here,  the  first  item  represents  the  useful  signal  and  the  second  one  with  the  spatial  and  temporal  £2 
frequencies  describes  the  amplitude  dynamic  distortions.  In  the  four-beam  geometry,  the  input  and  output  planes 
of  the  mitigation  system  always  coincide  whereas  in  the  two-beam  geometry,  they  are  separated  by  PRC  length. 
This  length  can  be  greater  than  the  diffraction  length.  That  is  why  we  recalculated  the  output  field  to  the  input 
plane.  In  a  real  experiment,  it  is  equivalent  to  addition  of  a  lens  to  form  the  input  image. 

In  parallel  with  (4),  in  four-beam  geometry,  we  characterized  the  quality  of  noise  mitigation  by  the  more 
convenient  value  -  the  overlapping  integral  H 

"Skri;|£j- 

2,3.  Theory 

To  illustrate  a  possibility  of  dynamic  noise  filtration  in  such  a  scheme,  we  can  use  the  simplest  case 
without  amplitude  distortions.  Let  as  signal  beam  with  the  amplitude  Ej  and  the  phase  distortions  (2.3)  as  plane 
pumping  wave  with  the  amplitude  A'?  are  directed  onto  PRC’s  input  plane 

+  .  (2  '7) 

Here  E  is  the  resulting  field  in  the  input  plane,  4  is  the  transverse  projection  of  the  wave  vector.  Under 
assumption  of  fixed  intensity  of  the  beams  and  with  ignoring  of  the  diffraction,  the  distribution  of  intensity  in  PRC 
can  be  found  as 


I=EE’=E]  +El+2E^E,J,{A^)f.yM,)exv\i[2k^  +mk)x]  + 

+2£|/?j  2  X-f„(^,y,(^,)exp{/|(2A,  +mk^+Ik}jx-nQl^. 

m=~<o  y^o 


(2.8) 


HqvqJJx)  is  the  m-ovdev  Bessel  function.  If  PRC  does  not  respond  to  the  high-frequency  components 
of  the  fringe  pattern  (2.8)  due  to  its  inertia,  the  gratings  at  the  spatial  frequencies  (Ik^+rnkJ  corresponding  to  the 
useful  information  will  be  only  written.  Thus,  in  the  case  of  so-called  diffusion  transportation  of  charge  ("nonlocal 
response",  Eo=0 ), 


oc  ^2  i{2k^  +  mk^)J ^^{A^)Qx^i{2k^  +mk^x 


(2.9) 


It  means  that  only  spectral  components  at  the  useful  spatial  frequencies  can  be  amplified.  When 
Q»  1  /  r  (r  is  the  PRC  characteristic  response  time)  and  spatial  frequencies  of  the  noise  and  the  useful  signal  are 
not  multiple,  one  can  obtain  steady-state  solution  as 


pou, 


(2.10) 


It  is  easy  to  see  that  depends  only  on  the  gain  constant  y.  It  means  that  the  larger  will  be  the 
difference  in  y  for  the  useful  and  noise  parts  of  A,  the  better  will  be  the  quality  of  filtration  of  dynamic 
distortions.  So,  in  two-beam  geometry,  the  maximal  magnitude  of  is  always  limited  by  the  maximal  value  of  y 
and  initial  noise  level. 

The  estimation  (2.10)  is  not  valid  for  the  four-beam  geometry.  Under  the  same  assumptions  and  for  the 
case  with  nonlocal  character  of  PRC  response, 

E°**'  oc  {Ik^  -  iAk^  cos(/r,x))exp|-/(A^x  -  A^  sin{A,x))| .  (211) 


Thus,  under  2k^  »  Ak^  that  is  typical  for  a  real  experiment,  the  output  beam  contains  only  the  useful 
information.  In  the  case  with  local  character  of  PRC  response,  the  useful  signal  can  be  completely  filtered  as  well. 
However,  even  in  these  simplest  examples,  the  solutions  change  radically  under  violation  of  conditions  mentioned 
above.  The  problem  becomes  much  more  complicated  if  we  try  to  take  into  account  diffraction,  self-action. 
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frequency  detuning,  real  PRC  thickness,  etc.  In  this  case,  it  is  impossible  to  find  both  steady-state  and  transient 
solutions  of  the  problem  in  analytical  form.  So,  computer  simulation  is  the  only  way  to  solve  this  problem. 

2.4.  Numerical  results 

All  parameters  of  our  simulation  corresponded  to  well-known  PRC  lnP;Fe  Its  electrooptical  constant 
is  1.45  pm/V,  /7=3.3,  and  the  free  path  length  of  carriers  is  determined  by  1.5*  10*'^  cm^/V.  We  represented  1 
mm  length  PRC  as  20-80  phase  screens  and  described  the  spatial  distributions  of  n  and  /  by  arrays  of  8192  points 
in  the  aperture  of  3  mm.  We  assumed  that  the  signal  beam  with  1  mm  diameter  formed  the  transmission  gratings 
with  5  pm  period. 

2.4. 1  ■  Two-beam  interaction 

For  this  geometry,  our  computer  simulation  confirmed  the  possibility  of  efficient  mitigation  of  phase 
dynamic  distortions.  Under  conditions  listed  in  section  1,  after  switching  the  mitigation  system  on,  decreases 
from  the  value  0''  up  to  the  steady-state  level  (2.9)  with  the  characteristic  time  r  (Fig.2.1b).  However,  with 
violation  of  any  conditions  listed  above,  the  dependence  changes  drastically.  Fig.2.1a  illustrates 

characteristic  oscillations  in  /f“^(?/r^for  the  case  when  ks=kn  and  Qx»  1.  These  oscillations  can  be  interpreted 
with  easy  by  simple  qualitative  reasoning.  In  this  case,  effective  energy  exchange  between  all  spectral  components 
of  the  signal  takes  place  because  the  refractive  index  gratings  contain  both  useful  and  noise  information.  That  is 
why  the  efective  filtration  of  dynamic  noise  does  not  occur.  The  range  of  k^,  where  we  can  observe  this  regime,  is 
determined  by  the  spectral  selectivity  of  written  3D  grating,  that  is,  by  the  aperture  of  interacting  beams. 

Nonlinear  response  of  PRC  is  usually  not  of  "pure”  nonlocal  character.  For  example,  such  a  character 
results  from  the  use  of  the  static  electric  field  Eq  to  increase  the  gain  constant  Our  simulation  has  shown  that 

strongly  depends  on  the  frequency  detuning  Aco  between  the  interacting  beams.  Fig. 2.2  represents  the 
steady-state  dependencies  for  different  relationships  between  the  local  and  nonlocal  parts  in  PRC 

nonlinear  response.  This  figure  illustrates  effects  arising  under  violation  of  the  Qi  <<  1  condition.  Because  our 
system  is  a  time  averaging  filter,  only  the  distortions  with  characteristic  times  that  are  less  than  r  can  be 
mitigated.  In  opposite  case,  the  gratings  at  frequency  £2  will  be  efficiently  written  in  PRC  that  leads  to  efficient 
amplification  of  the  noise  components  and  decreases  the  mitigation  quality.  Typical  dependencies  for 

some  values  of  the  parameter  Acol  are  shown  in  fig.2.3.  Fig. 2. 3a  corresponds  to  the  case  with  domination  of 
nonlocal  response.  Fig. 2. 3b  relates  to  the  opposite  case  when  the  quality  of  noise  mitigation  is  optimal  when 
Aru  =  0.  The  reason  of  this  result  is  following.  When  Aco  =0,  the  optimal  phase  shift  between  the  interference 
fringe  pattern  and  the  refractive  index  grating  {iifl)  is  achieved  for  some  shifted  spectral  components  of  the  noise. 
Therefore,  the  gain  constant  for  these  components  is  larger  than  the  gain  constant  for  the  useful  signal.  Obviously 
this  fact  will  decrease  the  mitigation  quality.  However,  one  can  optimize  the  phase  shift  by  means  of  some  well- 
known  ways.  It  is  possible  to  optimize  the  frequency  shift  of  pumping  in  relation  of  the  useful  signal  (see 
fig.. 2. 3b)  or  the  frequency  of  external  alternate  electric  field  The  last  technique  enables  one  to  realize  the 
highest  value  of  the  gain  constant. 


Fig.2.1.  versus  i/z  in  two-bcam  geometry:  =130 

cm'*,  kn  =130  (a)  and  100  (b)  enf*;  /1,=0.2; 
^;j=0.5;  Eo=0’,  Qr»\. 


Fig. 2. 2.  versus  T2rin  two-bcam  geometry:  k^^  =100 

cm'*;  k„=\3^  cm'*; /ls=0.2;  ^„=0.5;  Eo~  5.5  (a,b) 
and  15  (c)  kV/cm;  Acot=0{E)  and  1  (b,c). 
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Fig.2.3.  versus  t/z  and  Acoz  in  two-beam  geometry:  ks  =100  cm'^  /:„=130  cm'^;  >1^-0. 2;  >l«-0.5;  Eor^  (b)  and  5.5 

kV/cm  (a);  i2r»l. 

2.4.2.  Four-beam  interaction 

For  this  geometry  we  have  considered  the  problem  of  mitigation  of  both  phase  and  amplitude  dynamic 
distortions  and  studied  the  possibility  of  reconstruction  of  amplitude  image. 

According  to  our  computer  simulation  of  four-beam  geometry,  after  switching  the  system  on,  begins 
to  decrease  from  its  initial  value  to  the  steady-state  level  at  t  »  r  (fig.2.4a).  The  overlapping  integral 
dependence  //(//z)  exhibits  improvement  of  the  output  beam  quality,  too  (fig.2.4b).  However,  in  contrast  to  the 
two-beam  geometry  (see  fig.2.1,  curve  b),  the  characteristic  time  of  this  transient  process  is  determined  by  the 
new  time  constant  1/Q  that  is  much  smaller  than  r.  The  reason  of  such  a  difference  is  very  simple.  The  dynamic 
noise  penetrates  to  the  output  of  four-beam  mitigation  system  only  if  this  noise  has  been  written  in  the  refractive 
index  grating.  Therefore,  the  grating  amplitude  Sn  determines  only  the  energy  efficiency  of  the  system,  whereas 
the  accuracy  of  this  grating  determines  the  mitigation  quality.  In  our  case,  the  correct  grating  forms  during  the 
time  interval  1/n  and  after  that  the  output  signal  power  continue  to  grow  with  the  characteristic  time  r 
(fig.2.4c).  We  have  found  that  the  value  is  rather  sensitive  to  spatial  frequency  of  the  input  signal  (fig.2.5).  It 
means  that  the  limiting  level  of  distortion  mitigation  is  determined  by  the  small-scale  self-focusing  process.  The 
case  ks=kf,  is  followed  by  characteristic  oscillations  in  (fig.2.6).  However,  does  not  depend  on  the 

phase  shift  of  grating  and  frequency  detuning  between  the  interacting  beams  (fig.2.7). 


Fig. 2. 4.  (a),  H  (b),  and  (c)  versus  t/x  for  the  case  with  mitigation  of  phase  distortions  in  four-beam  geometiy:  ks 

=  100  cm'\  kr,=  m  cm'^Xv=0.2,  ^„=0.5;  Eo=0;  /2z»l. 


Our  simulation  has  confirmed  the  possibility  of  effective  mitigation  of  amplitude  dynamic  distortions  in 
four-beam  geometry  as  well  (fig.2.8  and  fig.2.9).  In  these  cases,  we  simulated  the  total  input  signal  as  sum  of  the 
amplitude  dynamic  noise  and  the  phase 


exp|-|^)  |(l  +  sin(A^x))  +  sin(A„x  - Q/) 

(2.12) 

or  amplitude 

^  o.5sin{A:^x))  +  5  sin(A  x  -  fl/) 

(2.13) 

71 


image,  respectively.  Hear,  a  is 
that  for  the  realizations  shown 


the  beam  radius,  £„  represents  the  dynamic  noise  amplitude.  It  is  vep-  important 
in  figures  8  and  9,  £„  is  five  times  larger  than  the  amplitude  of  the  useful  signal. 


Fig.2.5. 


/?“'  versus  kn  for  the  case  of  mitigation  of  Fig.2.6. 

r  St 

phase  distortions  in  four-beam  geometry: 

=130  cm-';  A,=0.2;  An=0.5;  Eo  =15  kV/cm; 

{2r»\, 


versus  t/r  for  Uic  case  of  miligalion  of  phase 
distortions  in  four-beam  geometry:  =  130  cm  , 
130  (a)  and  150  (b)  cm‘^;  4=0,2,  .4n=0.5;  Eo 
=  0;  Dr»l. 


Fig.2.7. 


/r'  versus  r/r  and  drur  for  the  case  of  mitigation  of  phase  distortions  in  four-beam  geometry:  h  -100  cm' 
£.=130  cm-';yl.=0.2;  rtn=0.5;  Eo=0  (b)  and  5.5  (a)  kV/cm;  /2r»l. 


Fig.2.9. 


H  tat  and  7°“'//"  (b)  versus  t/x  for  the  case  with  mitigation  of  amplitude  distortions  in  four-beam  geome^.  t/r  0 
corresponds  to  the  “switching  on”  moment;  £  =1 10  cm'',  k„=140  cm';  £3  =  0. 1;  £„  =  0.2, 0.6,  and  1;  £0  =  0;  Ox 
»1. 
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So,  our  simulation  has  proved  that  under  the  same  conditions,  in  comparison  with  the  systems  based  on 
the  two-beam  geometry,  the  mitigation  systems  based  on  the  four-beam  geometry  are  much  more  effective. 

Z5,  Conclusions  r*  u 

1.  A  new  promising  numerical  scheme  for  simulation  of  self-consistent  problems,  arising  from  the  problem  ot 

nonlinear  interaction  of  spatial-modulated  beams  in  PRC,  was  proposed.  ^  ^  ^ 

2.  The  scheme  was  tested  and  it  was  shown  that  two-  and  four-beam  geometry  of  nonlinear  interaction  in  PRC 

can  be  used  to  mitigate  phase  and  amplitude  dynamic  distortions.  •  •  dt?  r 

3.  The  scheme  enables  one  to  simulate  various  problems  of  dynamics  of  two-  and  four-beam  interaction  m  PRC 
with  taking  into  account  diffraction,  self-action,  exhaustion  of  pumping,  etc. 

4.  The  same  approach  can  be  extended  to  much  more  complicated  problems  of  nonlinear  interaction  such  as 

dynamics  of  self-pumping  PC-systems 
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